{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><font color=\"darkblue\">Logistic Regression</font></h2>\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Sigmoid function\n",
    "> $ \\displaystyle S(x) = \\frac{1}{1+e^{-x}} = \\frac{e^x}{1+e^x} \\qquad 0 < S(x) < 1 $\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/8/88/Logistic-curve.svg\" width=300 align=center>\n",
    "\n",
    "<p style=\"text-align:center\">(Picture from https://en.wikipedia.org/wiki/Sigmoid_function)</p>\n",
    "\n",
    "- Logit\n",
    "> $ \\displaystyle \\text{logit}(p) = \\log \\left( \\frac{p}{1-p} \\right) \\qquad 0 < p < 1 $\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/c/c8/Logit.svg\" width=300 align=center>\n",
    "\n",
    "<p style=\"text-align:center\">(Picture from https://en.wikipedia.org/wiki/Logit)</p>\n",
    "\n",
    "- What is the relationship between them?\n",
    "> $ \\displaystyle \\text{logit}(p) = \\log \\left( \\frac{p}{1-p} \\right) \\quad \\Longrightarrow \\quad e^{\\text{logit}(p)} = \\frac{p}{1-p} \\quad \\Longrightarrow \\quad p = \\frac{e^{\\text{logit}(p)}}{1+e^{\\text{logit}(p)}} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recall Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- General Form\n",
    "> $ \\displaystyle Y = X \\beta + \\epsilon $\n",
    "\n",
    "- Fit (least square)\n",
    "> $ \\displaystyle \\widehat{\\beta} = (X^\\top X)^{-1} X^\\top Y $\n",
    "\n",
    "- Predict\n",
    "> $ \\displaystyle \\widehat{Y} = X \\widehat{\\beta} $\n",
    "\n",
    "- Note that $ Y $ here is usually continuous. What if we have $ Y $ that is categorical, for example $ \\displaystyle Y \\in \\{ 0, 1 \\} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Features:  (100, 1)\n",
      "Target:  (100,)\n"
     ]
    }
   ],
   "source": [
    "# Generate data\n",
    "n_samples = 100\n",
    "np.random.seed(2018)\n",
    "\n",
    "x = np.random.normal(size=n_samples)\n",
    "y = (x > 0).astype(np.float)\n",
    "x[x > 0] *= 4\n",
    "x += .3 * np.random.normal(size=n_samples)\n",
    "x = x[:, np.newaxis]\n",
    "\n",
    "print('Features: ', x.shape)\n",
    "print('Target: ', y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAESNJREFUeJzt3XGMHOddxvHnuTtfnfO5Kdo9JPD5do3kAiYUJV5FCUFQyUFy2sr+A4FsKUiFqFbOTQkQAQ5FEQqyBBRVBWFKrRIqdY+kJhRkBQe3hCAEIsHnpk1jG1dXl9jnBuUa2oCoimvlxx+7d9k77+7M3u55fG++H2mknZl33/c3s+tHc+94dx0RAgCkZajoAgAAg0e4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABI0UtTA5XI5qtVqUcMDwLp0+vTpb0TERFa7wsK9Wq1qdna2qOEBYF2y/XKedkzLAECCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQoMxwt/2Y7Vdtv9Rhv23/ke052y/avm3wZQIAepHnyv1TknZ32X+PpO3N5YCkj/dfFtqZmZlRtVqVbY2MjMi2yuWyyuWyhoaGVK1WdfDgQVWr1aX1mZmZtn0dPHhQQ0NDsi3b2rx5s2ZmZnKN0a3PxeesZrnppps6Hku5XNbmzZuX2o6Pj2t8fPyaPkZGRnTw4MFl5yur7qy2MzMzKpfLS2OUy+WO56rTOL207VVW7XnGzTpXrfvzvh+ut15e7yJc9/oiInORVJX0Uod9n5C0v2X9vKTvy+pz586dgfzq9XqMjY2FpJ6WsbGxqNfry/qanp5u23ZoaCg2bNgw0D6LWnbt2nXN+WpXd6dzu9i2Xq/H6OjoNf0PDw93PFcrx+n22nWqqZ/3RWvtecbt1kee916/xzAIWcdQtEHWJ2k28uR2rkbdw/0pST/Rsv6MpFpWn4R7byqVyqqDrlKpLOtreHi47/Bciz6vx7Ky7m7ntlKprPq8t46T1Ue7mvp9X+SpfXHcbn3kfe/1cwyDkHUMRRtkfcoZ7m607c52VdJTEXFLm31PSfrdiPjn5vozkn4jIq754hjbB9SYutHU1NTOl1/O9RUJkDQ0NKQ8r1U7tvXGG28sW+/XWvR5PaysW+p8bhePaTXnvXWcrNeuXU159VP74rjd+ui2v13bomQdQ9EGWZ/t0xFRyxyzp17buyxpa8v6ZHPbNSLiaETUIqI2MZH5pWZoMTU1NbDnDg8P91vOmvR5PbQ7j53O7dTU1KrPe+vzsvoY5Gvbuj3vuN36yFtfP8cwCFnHULRC6stzea/u0zLvlfS0JEu6Q9K/5emTaZneMOfe28KcO3PuN0Jdi27IOXdJj0t6RdJ3Jc1Luk/S/ZLub+63pCOSvirpy8ox3x6E+6rU6/WlubvFOe5SqRSlUilsR6VSienp6ahUKkvrnd4809PTYXvpjTY+Pr4UCFljdOuzn7n3jRs3djyWUqkU4+PjS203bdoUmzZtahu609PTy85XVt1Zbev1epRKpaUxSqVSx3PVaZxe2vYqq/Y842adq9b9ed8P11svr3cRBlVf3nDPNee+Fmq1WvB97gDQm+s55w4AuMEQ7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACcoV7rZ32z5ve872oTb7p2w/a/sF2y/afs/gSwUA5JUZ7raHJR2RdI+kHZL2296xotlvSToWEbdK2ifpTwZdKAAgvzxX7rdLmouICxFxRdITkvauaBOS3t58fLOkrw+uRABAr/KE+xZJl1rW55vbWv22pHttz0s6IelD7TqyfcD2rO3ZhYWFVZQLAMhjUDdU90v6VERMSnqPpE/bvqbviDgaEbWIqE1MTAxoaADASnnC/bKkrS3rk81tre6TdEySIuJfJW2UVB5EgQCA3uUJ91OSttveZntUjRumx1e0uShplyTZ/mE1wp15FwAoSGa4R8RVSQ9IOinpnBr/K+aM7Udt72k2e0jSB2x/SdLjkt4fEbFWRQMAuhvJ0ygiTqhxo7R12yMtj89KumuwpQEAVotPqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJChXuNvebfu87Tnbhzq0+TnbZ22fsf0Xgy0TANCLkawGtoclHZH005LmJZ2yfTwizra02S7pYUl3RcQ3bX/vWhUMAMiW58r9dklzEXEhIq5IekLS3hVtPiDpSER8U5Ii4tXBlgkA6EWecN8i6VLL+nxzW6t3Snqn7X+x/Zzt3YMqEADQu8xpmR762S7p3ZImJf2T7R+NiG+1NrJ9QNIBSZqamhrQ0ACAlfJcuV+WtLVlfbK5rdW8pOMR8d2I+Jqkr6gR9stExNGIqEVEbWJiYrU1AwAy5An3U5K2295me1TSPknHV7T5GzWu2mW7rMY0zYUB1gkA6EFmuEfEVUkPSDop6ZykYxFxxvajtvc0m52U9Jrts5KelfRrEfHaWhUNAOjOEVHIwLVaLWZnZwsZGwDWK9unI6KW1Y5PqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AE5Qp327ttn7c9Z/tQl3Y/Yzts1wZXIgCgV5nhbntY0hFJ90jaIWm/7R1t2m2W9KCk5wddJACgN3mu3G+XNBcRFyLiiqQnJO1t0+53JP2epO8MsD4AwCrkCfctki61rM83ty2xfZukrRHxt906sn3A9qzt2YWFhZ6LBQDk0/cNVdtDkj4q6aGsthFxNCJqEVGbmJjod2gAQAd5wv2ypK0t65PNbYs2S7pF0j/a/g9Jd0g6zk1VAChOnnA/JWm77W22RyXtk3R8cWdEvB4R5YioRkRV0nOS9kTE7JpUDADIlBnuEXFV0gOSTko6J+lYRJyx/ajtPWtdIACgdyN5GkXECUknVmx7pEPbd/dfFgCgH3xCFQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQbnC3fZu2+dtz9k+1Gb/r9o+a/tF28/Yrgy+VABAXpnhbntY0hFJ90jaIWm/7R0rmr0gqRYR75L0pKTfH3ShAID88ly53y5pLiIuRMQVSU9I2tvaICKejYhvN1efkzQ52DIBAL3IE+5bJF1qWZ9vbuvkPklP91MUAKA/I4PszPa9kmqSfqrD/gOSDkjS1NTUIIcGALTIc+V+WdLWlvXJ5rZlbN8t6cOS9kTE/7XrKCKORkQtImoTExOrqRcAkEOecD8labvtbbZHJe2TdLy1ge1bJX1CjWB/dfBlAgB6kRnuEXFV0gOSTko6J+lYRJyx/ajtPc1mH5E0LukvbX/R9vEO3QEAroNcc+4RcULSiRXbHml5fPeA6wIA9IFPqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AE5Qp327ttn7c9Z/tQm/1vs/2Z5v7nbVcHXeiimZkZVatVDQ0NqVqtamZmJrN9uVyWbdlWuVy+5jlZfbbuL5fLS/0NDQ0t9dttydtuLftpd9wAEhYRXRdJw5K+KukHJI1K+pKkHSvaHJT0p83H+yR9JqvfnTt3Rq/q9XqMjY2FpKVlbGws6vV6x/ajo6PL2kuKDRs2LD0nq892+9frMjo62vFcAVgfJM1GRr5GRK5wv1PSyZb1hyU9vKLNSUl3Nh+PSPqGJHfrdzXhXqlU2oZWpVLpqX3rc7L67NbHelw6nSsA60PecM8zLbNF0qWW9fnmtrZtIuKqpNcllVZ2ZPuA7VnbswsLCzmGXu7ixYsD2d66L+u53fpYj1I7HgDtXdcbqhFxNCJqEVGbmJjo+flTU1MD2d66L+u53fpYj1I7HgDt5Qn3y5K2tqxPNre1bWN7RNLNkl4bRIGtDh8+rLGxsWXbxsbGdPjw4Y7tR0dHr9m+YcOGpedk9dlu/3o1Ojra8VwBSEzWvI0ac+gXJG3TmzdUf2RFmw9q+Q3VY1n9rmbOPaJxg7NSqYTtqFQqmTcI6/V6lEqlpTnnUql0zXOy+mzdXyqVlvqznWueO2+7teyn3XEDWH+Uc87djbbd2X6PpI+p8T9nHouIw7YfbQ5y3PZGSZ+WdKuk/5K0LyIudOuzVqvF7Oxs5tgAgDfZPh0Rtax2I3k6i4gTkk6s2PZIy+PvSPrZXosEAKwNPqEKAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCcn2IaU0GthckvVzI4O2V1fg2yxRwLDcmjuXGtN6OpRIRmV/OVVi432hsz+b51Nd6wLHcmDiWG1NKx9KKaRkASBDhDgAJItzfdLToAgaIY7kxcSw3ppSOZQlz7gCQIK7cASBBhHsL2x+x/e+2X7T917bfUXRNvbC92/Z523O2DxVdz2rZ3mr7WdtnbZ+x/WDRNfXL9rDtF2w/VXQt/bD9DttPNv+dnLN9Z9E1rZbtX2m+v16y/XjzdymSQbgv93lJt0TEuyR9RdLDBdeTm+1hSUck3SNph6T9tncUW9WqXZX0UETskHSHpA+u42NZ9KCkc0UXMQB/KOnvIuKHJP2Y1ukx2d4i6Zck1SLiFjV+iGhfsVUNFuHeIiI+FxFXm6vPqfF7sevF7ZLmIuJCRFyR9ISkvQXXtCoR8UpEfKH5+H/UCJAtxVa1erYnJb1X0ieLrqUftm+W9JOS/kySIuJKRHyr2Kr6MiLppubvPo9J+nrB9QwU4d7ZL0p6uugierBF0qWW9Xmt40BcZLuqxs83Pl9sJX35mKRfl/RG0YX0aZukBUl/3pxi+qTtTUUXtRoRcVnSH0i6KOkVSa9HxOeKrWqw3nLhbvvvm3NsK5e9LW0+rMbUwExxlcL2uKS/kvTLEfHfRdezGrbfJ+nViDhddC0DMCLpNkkfj4hbJf2vpHV5b8f296jxl+02Sd8vaZPte4utarBy/YZqSiLi7m77bb9f0vsk7Yr19f9EL0va2rI+2dy2LtneoEawz0TEZ4uupw93SdrT/JH5jZLebrseEesxSOYlzUfE4l9RT2qdhrukuyV9LSIWJMn2ZyX9uKR6oVUN0Fvuyr0b27vV+PN5T0R8u+h6enRK0nbb22yPqnFz6HjBNa2Kbasxr3suIj5adD39iIiHI2IyIqpqvCb/sE6DXRHxn5Iu2f7B5qZdks4WWFI/Lkq6w/ZY8/22S+v05nAnb7kr9wx/LOltkj7feL31XETcX2xJ+UTEVdsPSDqpxp3/xyLiTMFlrdZdkn5e0pdtf7G57Tcj4kSBNaHhQ5JmmhcQFyT9QsH1rEpEPG/7SUlfUGMK9gUl9klVPqEKAAliWgYAEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQoP8HD0dAzrFvbm4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y, color='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For this example, if we fit a linear regression model\n",
    "> $ \\displaystyle y = \\beta_0 + \\beta_1 x + \\epsilon $\n",
    ">\n",
    "> $ \\displaystyle \\widehat{y} = \\widehat{\\beta}_0 + \\widehat{\\beta}_1 x $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit a linear regression model\n",
    "XX = np.ones((x.shape[0], 2))\n",
    "XX[:,1] = x.ravel()\n",
    "bHat = linalg.inv(XX.T.dot(XX)).dot(XX.T).dot(y)\n",
    "b0, b1 = bHat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4lOW5x/HvnYQAAWRJ2CEJCC4sbkRAhGjrUcG1i+1lXaqtFgUX4CyCyzm1Vs+xrUewLrW41GpibSWgFNlUVLDWJbgkLLLvoGyyyBaSuc8fGTghJCQkk7yTmd/nuuZi3pln3ueemfDLk3dm7jF3R0RE4ktC0AWIiEj9U/iLiMQhhb+ISBxS+IuIxCGFv4hIHFL4i4jEIYW/iEgcUviLiMQhhb+ISBxKCrqAyqSlpXlmZmbQZYiINCjz58/f6u5tqxoXteGfmZlJfn5+0GWIiDQoZramOuN02EdEJA4p/EVE4pDCX0QkDin8RUTikMJfRCQOKfxFROJQrcPfzLqa2TtmtsjMFprZqArGmJn93syWm1mBmZ1V23lFRKTmIrHyLwb+zd17AQOB28ysV7kxw4Ce4dNw4A8RmFdEJOZs2X2ANdv21Pk8tQ5/d9/k7p+Gz+8GFgOdyw27EnjRS30ItDKzjrWdW0QkVrg7kz9dz7DH5vL+8q11Pl9EP+FrZpnAmcBH5a7qDKwrs70+fNmmSM4vItIQbdixj3unFPL1rgP86cb+9O3Sss7njFj4m1lzIA8Y7e67ariP4ZQeFiI9PT1SpYmIRKVQyMn9aA3j31rGz8/N5JbzTqRRYv28Dyci4W9mjSgN/lx3n1zBkA1A1zLbXcKXHcHdJwITAbKysjwStYmIRKOVW75lXF4hxaEQf7tlID3atajX+Wsd/mZmwHPAYnd/tJJhU4HbzewVYACw0911yEdE4k5xSYhn5q1i4twV3HlBT356TiaJCVbvdURi5X8ucD1QaGafhy+7B0gHcPengenAJcByYC/wswjMKyLSoCzcuJOxeQW0Tklm6u2D6domJbBaah3+7v4+cMxfW+7uwG21nUtEpCHaf7CEx+cs45WP1zFu2Clc1a8LpQdNghO1/fxFRGLB/DXbuWtSAT3aNWfGqCG0O6FJ0CUBCn8RkTqx50Axv5u1hOmFm/jVFb0Z1je6Ptqk8BcRibC5S7dwz5RCBnRLZfaYbFqlJAdd0lEU/iIiEbJz70F+/cYi/rliG//9g76cd1KVX6UbGHX1FBGJgJkLNnHRhPdolpzIrDHZUR38oJW/iEitbN69n1++vpAlX+/miWvO4uzMNkGXVC1a+YuI1IC782r+OoZNmEe3tGZMv3NIgwl+0MpfROS4rdu+l3umFLLt2yL+/PP+9Olc943YIk3hLyJSTaGQ8+I/V/PY28u4eUh3hmd3r7dGbJGm8BcRqYblm79lXF4BAK/eOoge7ZoHXFHtKPxFRI7hYEmIiXNX8uy8lYy58CSuG5BBQgCN2CJN4S8iUokFG3Zy16QC0lo05u93DKZL6+AasUWawl9EpJz9B0t47O1lvJq/jruHncoPzuoceCO2SFP4i4iU8cnq7YydVMApHVswY1Q2bVs0DrqkOqHwFxEBvj1QzG9nfsmshV/xqyt6M7RPdDViizSFv4jEvXeXbObeKQsYdGIqs0efR8uURkGXVOcU/iISt77ZU8Sv31jEx6u28/AP+zKkZ3T344mkhvnpBBGRWnB33ijYxEUT5tKyaSNmjc6Oq+AHrfxFJM5s3rWf+15bwMqte3j6un70y2gddEmB0MpfROKCu/O3T9Yx7LF5nNyhBW/cOThugx+08heROLBu+17unlzIzn0HeemmAfTqdELQJQVO4S8iMask5Pz5g9U8PmcZt5x3IjcP7kZSA23EFmkKfxGJScu+3s1deQU0Skwgb8Qgurdt2I3YIk3hLyIxpag4xNPvreCFD1bzrxeexDX902OiEVukKfxFJGYUrN/BXZMK6NCyCdPuGEynVk2DLilqKfxFpMHbf7CE8W8uJe/T9dx3aS+uPKNTzDViizSFv4g0aB+u3Ma4vAL6dG7JzNHZpDWPzUZskabwF5EGaff+gzw840veXryZB67szUW9OwRdUoOi9zyJSIMz58uvuXj8XELuzBqTreCvAa38RaTB2L6niAf+vpBP1+7gkR+dzqAeaUGX1GBp5S8iUc/dmfrFRi4aP5e05o2ZOXqIgr+WIrLyN7PngcuAze7ep4LrzwdeB1aFL5rs7g9EYm4RiW1f7SxtxLZm2x6e+Wk/zkyP3348kRSpwz4vAE8ALx5jzDx3vyxC84lIjHN3XvlkHb+btYTrB2bw5LVn0jgpMeiyYkZEwt/d55pZZiT2JSKyZtsexuUVsreomJd/MYBTOqgRW6TV5wu+55jZF8BG4N/dfWH5AWY2HBgOkJ6eXo+liUg0KAk5f/rHKp58Zzkjz+/Bzwd3I1GtGepEfYX/p0CGu39rZpcArwE9yw9y94nARICsrCyvp9pEJAos+aq0EVvTRglMGXkumWnNgi4pptVL+Lv7rjLnp5vZU2aW5u5b62N+EYleRcUhnnp3OS/+cw3/ftHJXH12VzViqwf1Ev5m1gH42t3dzPpT+hbTbfUxt4hEr8/X7WDspAK6tG7KG3cOpmNLNWKrL5F6q+dfgPOBNDNbD/wSaATg7k8DVwEjzKwY2Adc7e46rCMSp/YVlfDom0uY8tlG/vOyU7nidDViq2+RerfPT6q4/glK3woqInHugxVbuXtyIWd0bcWs0UNIVSO2QKi9g4jUi137D/I/0xfz7pItPPi9PlxwavugS4prau8gInXurUWljdjMjFljshX8UUArfxGpM9u+PcD9f19EwfodPPrjMzjnxNSgS5IwrfxFJOLcndc/38DFE+bRsWUTZo7KVvBHGa38RSSiNu7Yx32vLWDjjn08d0MWp3dtFXRJUgGFv4hERCjkvPzxWh59cyk3Dsrk6ev6kZykgwvRSuEvIrW2ausexuUVcKA4xCvDB3JS+xZBlyRVUPiLSI0Vl4R47v1VPP3eCm7/bk9uHJSpRmwNhMJfRGpk8aZdjM0roEWTJF6/bTDpqSlBlyTHQeEvIsflQHEJT85ZTs5Haxk79GR+nNVVrRkaIIW/iFTbp2u/YeykAjLTmjFj1BDan9Ak6JKkhhT+IlKlvUXFPDJrKX8v2MgvL+/FpX07arXfwCn8ReSY3l+2lbunFHB2Rhtmj86mdbPkoEuSCFD4i0iFdu47yENvLOIfy7fx4Pf78J2T2wVdkkSQPoEhIkeZtfArLhr/Ho2TEpk5eoiCPwZp5S8ih23ZfYD7py5k0aZd/P7qMxnQXf14YpVW/iKCuzP50/UMe2wuXdukMGPUEAV/jNPKXyTObdixj3smF7J59wH+dGN/+nZpGXRJUg8U/iJxKhRycj5aw/g3l3LzkO4Mz+5Oo0QdDIgXCn+ROLRiy7eMyyugJOS8eus59GinRmzxRuEvEkeKS0JMnLeSZ+auZNQFPbn+HDVii1cKf5E4sXDjTsbmFdA6JZmptw+maxs1YotnCn+RGLf/YAmPz1nGKx+vY9ywU7iqXxe1ZhCFv0gsy1+9nbvyCjipXQtmjBpCOzVikzCFv0gM2nOgmN/NWsL0wk386oreDOvbMeiSJMoo/EVizNylW7h7ciEDu6cye0w2rVLUiE2OpvAXiRE79hbx4BuL+eeKbfz3D/py3kltgy5Jopg+0SESA2YUbuKi8XNplpzIrDHZCn6pklb+Ig3Y5t37+eXrC1ny9W6evPYszs5sE3RJ0kBo5S/SALk7r+avY9iEeXRv24zpdw5R8MtxicjK38yeBy4DNrt7nwquN+Ax4BJgL3Cju38aiblF4s267Xu5Z0oh2/cU8eJN/endSY3Y5PhFauX/AjD0GNcPA3qGT8OBP0RoXokiubm5ZGZmYmYkJCRgZpgZaWlpjBw5kszMTBISEsjMzCQ3N/eI244cOfKI27Ro0YKRI0eSlpZ2+LLmzZuTlpZW4T5yc3OPGFvdU6NGjQ7vMy0t7fA+EhMTj3m7tLQ0cnNzD9/nyu5XZdeXr/fQY1T+srL7eyknl8wLb2DQ/a8x7blHmT7mO1w66PSj5qzu81Sdmso/xpXd17LPfVJSEmZW4eNRV6p6HmJ9/hpx94icgExgQSXX/RH4SZntJUDHY+2vX79+Lg1HTk6Op6SkOFCtU0pKiufk5Li7+4gRI6p9u4r2kZOT440aNarRPmpzSkhI8OTk5ErvV0WPSUpKio8YMeKo21V2Sk5O9pycHH/02Ze90/WPePtrfuNJbTpXOmdNnqdDNVX0GB6av7LbVXZdTWqrqWPVVh+Cnr88IN+rk9nVGVStHR07/KcBg8tsvw1kHWt/Cv+GJSMj47jDMyMjw93dExMTaxzAGRkZNZq7Lk+H7ldldR3X/U1I9Ixhwz1j9Cve/MxLHeyYc9b0eTpWTcd6jKvz+Fe3tpo6Vm31Iej5y6tu+Fvp2Nozs0xgmld8zH8a8LC7vx/efhsY6+755cYNp/SwEOnp6f3WrFkTkdqk7iUkJHC8P0tmRigUqlWfmUO3jdTPcSQcul81eUzKSm5/IqnDRlGy5xu2z36S4p2bq5yzKjV9nqDix7g6j391a6upyu5TXc8bLfNXMO98d8+qalx9vdtnA9C1zHaX8GVHcPeJ7p7l7llt2+p9yg1Jenp6jW+TmJhYq3lrMnddOlRPZXVVdX8tKZlW591Aux/dz65PptD04+fp3Kppteasbm3HU9OxHuPqPP51/fwcq7b6EPT8NVadPw+qc+LYh30uBWYABgwEPq5qfzrs07DomP/R96smx/wbd+ntnW5+2tOuHOsJKa2Oecy9ojlr8jzpmH/Dnr886vOYP/AXYBNwEFgP3ATcCtwavt6AJ4EVQCFVHO93hX+DlJOTc/j4p9n/H5tOTU31ESNGeEZGhpuZZ2RkHPUfY8SIEUfcpnnz5j5ixAhPTU09fFmzZs08NTW1wn3k5OQcMba6p6SkpMP7TE1NPbyPhISEY94uNTX1cPAd635Vdn35elM7dPYhY57wrre/5E17nnPEHBU9voeO0Vc0Z3WfpyprqmT+yh7/SNRWU1U9D7E+f1nVDf+IHfOPtKysLM/Pz696oEgD986Szdw3ZQHn9kjl3kt60TKlUdAlSQNW3WP+au8gEpBv9hTx62mL+Hj1dn7zw9MY3DMt6JIkjqi9g0g9c3feKNjExRPm0iolmVmjsxX8Uu+08hepR1/v2s9/vraAlVv38Ifr+tEvo3XQJUmcUviL1AN352/56/jtzCVcOyCdx685k8ZJNX+Lq0htKfxF6tjabXu5e0oBu/YV89JNA+jV6YSgSxJR+IvUlZKQ88IHq3lizjJuPe9EbhrcjaREvcwm0UHhL1IHln29m7vyCmiUmEDeiEF0b9s86JJEjqDwF4mgouIQT7+3ghc+WM2/XngS1/RPJyGh5r2LROqKwl8kQr5Yt4OxeQV0bNmEaXcMplMV/XhEgqTwF6mlfUUlTHhrKXmfrue+S3tx5RmdatWpVKQ+KPxFauHDldsYl1dA3y6tmDk6m7TmjYMuSaRaFP4iNbB7/0EenvElby/ezK+/14cLe7UPuiSR46L3nYkcpzlffs3F4+cScmfWmGwFvzRIWvmLVNO2bw/wwLRFfLZ2B4/86HQG9VA/Hmm4tPIXqYK7M/WLjVw8YR5tmzdm1uhsBb80eFr5ixzDVzv3c99rhazdvpdnftqPM9PViE1ig8JfpALuziufrON3s5Zw/cAMnrq2H8lJ+kNZYofCX6ScNdv2MC6vkL1Fxbz8iwGc0kGN2CT2KPxFwkpCzp/+sYon31nOyPN78PPB3UhUawaJUQp/EWDJV6WN2Jo2SmDKyHPJTGsWdEkidUrhL3GtqDjEk+8s56UP1/AfF5/M1Wd3VWsGiQsKf4lbn6/bwV2TviC9TQrT7xxCh5ZNgi5JpN4o/CXu7Csq4X9nL+G1zzfyX5f34vLTOmq1L3FH4S9x5YMVWxmXV8iZ6a2YPSabNs2Sgy5JJBAKf4kLO/cd5OEZi3l3yRYe/F4fLjhV/XgkvulTKxLz3lxU2ogtwYxZY7IV/CJo5S8xbOu3B7h/6kIWbNjJhKvPYGD31KBLEokaWvlLzHF3XvtsA0MnzKVzq6bMGJWt4BcpRyt/iSkbd+zjvtcWsHHHPp6/8WxO69Iq6JJEopLCX2JCKOS8/PFaHn1zKTcOyuTp69SITeRYFP7S4K3auoexeQUUFYd4ZfhATmrfIuiSRKJeRJZGZjbUzJaY2XIzG1fB9Tea2RYz+zx8ujkS80p8Ky4J8fR7K/jBU//g4t4dyBsxSMEvUk21XvmbWSLwJHAhsB74xMymuvuickP/6u6313Y+EYBFG3cxNq+AE5om8fptg0lPTQm6JJEGJRKHffoDy919JYCZvQJcCZQPf5FaO1BcwhNzlvPyR2sZO/QUfpTVRa0ZRGogEuHfGVhXZns9MKCCcT80s2xgKTDG3ddVMEakUvPXfMPYvAK6pTVj+qghtD9BjdhEaqq+XvD9O/AXdz9gZrcAfwa+W36QmQ0HhgOkp6fXU2kS7fYWFfO7WUuYVrCJ+y/vzSV9O2i1L1JLkXjBdwPQtcx2l/Blh7n7Nnc/EN58FuhX0Y7cfaK7Z7l7Vtu2bSNQmjR07y/bysUT5rJz70Fmj87mUnXgFImISKz8PwF6mlk3SkP/auCasgPMrKO7bwpvXgEsjsC8EsN27j3IQ9MX8Y/l23jw+334zsntgi5JJKbUOvzdvdjMbgdmAYnA8+6+0MweAPLdfSpwp5ldARQD24EbazuvxK6ZC77il1MXcHHvDswak03zxvo4ikikmbsHXUOFsrKyPD8/P+gypB5t2V3aiG3Rpl385oen0b9bm6BLEmlwzGy+u2dVNU6ff5fAuTt589cz7LG5pKemMGPUEAW/SB3T39MSqA079nHP5EI27z7An27sT98uLYMuSSQuKPwlEKGQk/PRGia8tYybBndjeHZ3GiXqD1GR+qLwl3q3Ysu3jMsrIOTwt1vOoUe75kGXJBJ3FP5Sbw6WhHhm3kqembuSURf05KfnZJKQoPfsiwRB4S/1YsGGnYzNK6BNs2Sm3j6Yrm3UiE0kSAp/qVP7D5bw+JxlvPLxOsYNO4Wr+qkRm0g0UPhLnclfvZ278go4qV0LZoweQrsWasQmEi0U/hJxew6UNmKbXriJX13Rm2F9OwZdkoiUo/CXiHpv6RbumVzIOSemMntMNq1SkoMuSUQqoPCXiNixt4hfT1vMhyu38T8/6Ev2SerKKhLN9KkaqbUZhZu4aPxcWjRJYvaYbAW/SAOglb/U2OZd+/mv1xeydPNunrr2LLIy1Y9HpKHQyl+Om7vzav46hj02jxPbNWP6nUMU/CINjFb+clzWbd/LPVMK2b6niBdv6k/vTmrEJtIQKfylWkpCzov/XM3v317G8OwT+cWQbiSpEZtIg6Xwlyot37ybsXmFJBhMGjGIE9uqEZtIQ6fwl0odLAnxx/dW8Nz7q/jXC0/i2gEZasQmEiMU/lKhBRt28h+TCmjXojF/v2MwXVqrEZtILFH4yxH2HyxhwlvLmDR/HfdccirfP7OzGrGJxCCFvxz28artjMsr4NROJzBjVDZtWzQOuiQRqSMKf2H3/oP8duYSZi/6igeu7MPFvTsEXZKI1DG9Vy/OvbNkM0MnzKOoOMTs0ecp+EXihFb+ceqbPUX8etoiPlmznd/88DQG90wLuiQRqUda+ccZd2dawUYumjCXVinJzBqdreAXiUNa+ceRr3ft5z9fW8DKrXt4+rp+9MtoHXRJIhIQhX8ccHf+lr+O385cwrUD0nn8mjNpnJQYdFkiEiCFf4xbu20v4yYXsHt/MTk3D+DUjicEXZKIRAGFf4wqCTkvfLCaJ+Ys49bzTuSmwWrEJiL/T+Efg5Z+vZu7JhWQnJTA5JHn0i2tWdAliUiUichS0MyGmtkSM1tuZuMquL6xmf01fP1HZpYZiXnlSEXFIX7/9jKunvghV/Xrwiu/GKjgF5EK1Tr8zSwReBIYBvQCfmJmvcoNuwn4xt17AOOB39R23liRm5tLZmYmCQkJZGZmkpubW+m4tLQ0zAwzIy0t7fDY3NxcMs/KJnP4Uzz658nc1GUL1w0s7cB5aP9mRmJi4uHbR8up7P0QkXrk7rU6AecAs8ps3w3cXW7MLOCc8PkkYCtgx9pvv379PNbl5OR4SkqKA4dPKSkpnpOTc9S45OTkI8YB3qhRIx8+4jZve+Fw73LbS55y6nlH7KOi/UfjKTk5+aj7LCI1A+R7NbLbSsfWnJldBQx195vD29cDA9z99jJjFoTHrA9vrwiP2VrZfrOysjw/P79WtUW7zMxM1qxZc9TlGRkZrF69uspxjbv2JXXYHRRtWsb2t/5IaN+uI/YBVHi7aFT+PotIzZjZfHfPqmpcVL3ga2bDgeEA6enpAVdT99auXVuty8tvW3IKrc//GU17nM322U+xb/nH1d53tGpo9Yo0dJF4wXcD0LXMdpfwZRWOMbMkoCWwrfyO3H2iu2e5e1bbtm0jUFp0q+wXXPnLy2437Z5Fp5ueBDM2PjuSolXzK91HQ/oF2pBqFYkFkQj/T4CeZtbNzJKBq4Gp5cZMBW4In78KmOO1Pd4UAx566CFSUo78hqyUlBQeeuiho8Y1aZlG2mX/Tut/uYVtb4xn+6wnSPKDDB8+vNJ9VLT/aJScnHzUfRaROladFwaqOgGXAEuBFcC94cseAK4In28CvAosBz4Gule1z3h4wde99MXcjIwMNzPPyMg46oXPUCjkr3223nvfO9U7DrvNLamxA56amnp47LH2ceg6wBMSEgJ/cbf8qez9EJHao75e8K0r8fCCb1U27dzHfVMWsP6bffzmqtM4o2uroEsSkSjXIF/wlVKhkPPKJ+t4ZPYSfnpOBn+4rh/JSWrNICKRo/CPMqu37mHc5AL2HQzxl18M5OQOLYIuSURikMI/SpSEnOffX8VT7y7ntu/04GfndiMxwYIuS0RilMI/Ciz5ajd3TfqClOQkXrvtXDJS1Y9HROqWwj9AB4pLeOqdFbz04Rr+4+KTufrsrphptS8idU/hH5DP1n7D2LwC0tukMP3OIXRo2STokkQkjij869neomL+d/ZSXv98I7+8vBeXndZRq30RqXcK/3r0wfKtjJtcyFnprZg9Jps2zZKDLklE4pTCvx7s3HeQ/5m+mLlLt/Dg9/vw3VPaB12SiMQ5fXKojr256GsuHj+XxARj1phsBb+IRAWt/OvI1m8PcP/UhSzYsJMJV5/BwO6pQZckInKYVv4R5u5M+Ww9QyfMpXPrpswcna3gF5Goo5V/BG3csY97pxSyaed+nr/xbE7rokZsIhKdFP4REAo5uR+vZfybS/nZoExuOe9ENWITkaim8K+lVVv3MDavgIMlIf46fCA926sRm4hEP4V/DRWXhHj2/VX88b0V3PHdntwwKFON2ESkwVD418CijbsYm1dAy6aNmHr7YLq2if6vShQRKUvhfxwOFJfwxJzlvPzRWsYOPYUfZXVRawYRaZAU/tU0f01pI7buac2YPmoI7U9QIzYRabgU/lXYc6CYR2YvYVrBJu6/vDeX9O2g1b6INHgK/2OYt2wLd08upH+3NswenU1rNWITkRih8K/Azr0HeWj6Iv6xfBsPfr8P3zm5XdAliYhElD6JVM7MBV9x0YT3aNIokVljshX8IhKTtPIP27x7P/dPXciXm3bz+E/Oon+3NkGXJCJSZ+J+5e/u5M1fzyWPzSMjtfSdPAp+EYl1cb3yX//NXu6ZsoCtuw/wws/606dzy6BLEhGpF3EZ/qGQk/PRGsa/uZSbh3RneHZ3GiXG/R9BIhJH4i78V2z5lnF5BYQcXr11ED3aNQ+6JBGRehc34X+wJMQz81by7LxVjLqgJ9cPzCBBjdhEJE7FRfgv2LCTsXkFtGmWzOu3natGbCIS92I6/PcfLOH3by/jr5+s4+5LTuWHZ3VWawYREWoZ/mbWBvgrkAmsBn7s7t9UMK4EKAxvrnX3K2ozb3Xkr97OXXkFnNy+BTNGD6FdCzViExE5pLYr/3HA2+7+sJmNC2+PrWDcPnc/o5ZzVYu788C0RbxRsIlfXdGbYX071se0IiINSm3D/0rg/PD5PwPvUnH41xsz46z01oy6oCetUtSITUSkIrV9c3t7d98UPv8V0L6ScU3MLN/MPjSz79VyzipdfnonBb+IyDFUufI3s7eADhVcdW/ZDXd3M/NKdpPh7hvMrDswx8wK3X1FBXMNB4YDpKenV1m8iIjUTJXh7+7/Utl1Zva1mXV0901m1hHYXMk+NoT/XWlm7wJnAkeFv7tPBCYCZGVlVfaLREREaqm2h32mAjeEz98AvF5+gJm1NrPG4fNpwLnAolrOKyIitVDb8H8YuNDMlgH/Et7GzLLM7NnwmFOBfDP7AngHeNjdFf4iIgGq1bt93H0bcEEFl+cDN4fPfwD0rc08IiISWWplKSIShxT+IiJxSOEvIhKHzD0631FpZluANUHXUUYasDXoIqoQ7TVGe30Q/TVGe30Q/TVGe31Quxoz3L1tVYOiNvyjjZnlu3tW0HUcS7TXGO31QfTXGO31QfTXGO31Qf3UqMM+IiJxSOEvIhKHFP7VNzHoAqoh2muM9vog+muM9vog+muM9vqgHmrUMX8RkTiklb+ISBxS+NeAmf2bmXm4UV3UMLPfmdmXZlZgZlPMrFXQNR1iZkPNbImZLQ9/61vUMLOuZvaOmS0ys4VmNiromipjZolm9pmZTQu6lvLMrJWZTQr/DC42s3OCrqk8MxsTfo4XmNlfzCzw73c1s+fNbLOZLShzWRsze9PMloX/bR3peRX+x8nMugIXAWuDrqUCbwJ93P00YClwd8D1AKWBBTwJDAN6AT8xs17BVnWEYuDf3L0XMBC4LcrqK2sUsDjoIirxGDDT3U8BTifK6jSzzsCdQJa79wESgauDrQqAF4Ch5S479BW5PYG3w9sRpfA/fuOBu4Coe7HE3We7e3F480OgS5D1lNEfWO7uK929CHiF0q8AjQruvslvHzXEAAAClElEQVTdPw2f301paHUOtqqjmVkX4FLg2arG1jczawlkA88BuHuRu+8ItqoKJQFNzSwJSAE2BlwP7j4X2F7u4isp/Wpcwv9G/BsQFf7HwcyuBDa4+xdB11INPwdmBF1EWGdgXZnt9URhuAKYWSalXzb0UbCVVGgCpQuPUNCFVKAbsAX4U/iw1LNm1izoosoKf6nUI5T+1b4J2Onus4OtqlLV/YrcGlP4l2Nmb4WPB5Y/XQncA/xXFNd3aMy9lB7KyA2u0obHzJoDecBod98VdD1lmdllwGZ3nx90LZVIAs4C/uDuZwJ7qINDFbURPm5+JaW/qDoBzczsumCrqpqXviUz4kcaatXPPxZV9rWVZtaX0h+aL8wMSg+pfGpm/d39q6DrO8TMbgQuAy7w6Hkf7waga5ntLuHLooaZNaI0+HPdfXLQ9VTgXOAKM7sEaAKcYGY57h4t4bUeWO/uh/5imkSUhT+lXzi1yt23AJjZZGAQkBNoVRWr1lfk1oZW/tXk7oXu3s7dM909k9If9rPqM/irYmZDKT0scIW77w26njI+AXqaWTczS6b0RbapAdd0mJX+Nn8OWOzujwZdT0Xc/W537xL+2bsamBNFwU/4/8E6Mzs5fNEFRN/Xta4FBppZSvg5v4Aoe1G6jCq/Ire2tPKPLU8AjYE3w3+dfOjutwZbErh7sZndDsyi9B0Wz7v7woDLKutc4Hqg0Mw+D192j7tPD7CmhugOIDf8C34l8LOA6zmCu39kZpOATyk9LPoZUfBpXzP7C3A+kGZm64FfUvqVuH8zs5so7W7844jPGz1HBkREpL7osI+ISBxS+IuIxCGFv4hIHFL4i4jEIYW/iEgcUviLiMQhhb+ISBxS+IuIxKH/A++Feny9y/buAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_test = np.linspace(-5, 10, 300)\n",
    "\n",
    "plt.plot(X_test, b0 + b1 * X_test, linewidth=1);\n",
    "plt.scatter(x, y, color='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression (Binary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- General Form\n",
    "> $ \\displaystyle \\text{logit} \\left[ \\pi(x) \\right] = \\log \\left[ \\frac{\\pi(x)}{1-\\pi(x)} \\right] = \\beta_0 + \\beta_1 x \\qquad $ where $ \\displaystyle \\qquad \\pi(x) = \\mathbb{P}(y=1 \\ | \\ X = x) $\n",
    "\n",
    "- Fit (maximum likelihood)\n",
    "> $ \\displaystyle \\mathcal{L}(\\beta_0, \\beta_1) = \\prod\\limits_{i: y_i = 1} \\pi(x_i; \\beta_0, \\beta_1) \\prod\\limits_{j: y_{j} = 0} [1 - \\pi(x_j; \\beta_0, \\beta_1)] $\n",
    ">\n",
    "> $ \\displaystyle l(\\beta_0, \\beta_1) = \\log \\mathcal{L}(\\beta_0, \\beta_1) = \\sum_{i=1}^{n} \\left[y_i \\log \\pi(x_i; \\beta_0, \\beta_1) + (1-y_i) \\log(1-\\pi(x_i; \\beta_0, \\beta_1)) \\right] = \\sum_{i=1}^{n} \\left[y_i (\\beta_0 + \\beta_1 x_i) - \\log (1 + \\exp (\\beta_0 + \\beta_1 x_i))  \\right] $\n",
    ">\n",
    "> $ \\displaystyle \\widehat{\\beta}_0, \\widehat{\\beta}_1 = \\arg \\max_{\\beta_0, \\beta_1} l(\\beta_0, \\beta_1) $\n",
    "\n",
    "- Predict\n",
    "> $ \\displaystyle \\widehat{\\pi}(x) = \\frac{\\exp \\left( \\widehat{\\beta}_0 + \\widehat{\\beta}_1 x \\right)}{1 + \\exp \\left( \\widehat{\\beta}_0 + \\widehat{\\beta}_1 x \\right) } $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAESNJREFUeJzt3XGMHOddxvHnuTtfnfO5Kdo9JPD5do3kAiYUJV5FCUFQyUFy2sr+A4FsKUiFqFbOTQkQAQ5FEQqyBBRVBWFKrRIqdY+kJhRkBQe3hCAEIsHnpk1jG1dXl9jnBuUa2oCoimvlxx+7d9k77+7M3u55fG++H2mknZl33/c3s+tHc+94dx0RAgCkZajoAgAAg0e4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABI0UtTA5XI5qtVqUcMDwLp0+vTpb0TERFa7wsK9Wq1qdna2qOEBYF2y/XKedkzLAECCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQoMxwt/2Y7Vdtv9Rhv23/ke052y/avm3wZQIAepHnyv1TknZ32X+PpO3N5YCkj/dfFtqZmZlRtVqVbY2MjMi2yuWyyuWyhoaGVK1WdfDgQVWr1aX1mZmZtn0dPHhQQ0NDsi3b2rx5s2ZmZnKN0a3PxeesZrnppps6Hku5XNbmzZuX2o6Pj2t8fPyaPkZGRnTw4MFl5yur7qy2MzMzKpfLS2OUy+WO56rTOL207VVW7XnGzTpXrfvzvh+ut15e7yJc9/oiInORVJX0Uod9n5C0v2X9vKTvy+pz586dgfzq9XqMjY2FpJ6WsbGxqNfry/qanp5u23ZoaCg2bNgw0D6LWnbt2nXN+WpXd6dzu9i2Xq/H6OjoNf0PDw93PFcrx+n22nWqqZ/3RWvtecbt1kee916/xzAIWcdQtEHWJ2k28uR2rkbdw/0pST/Rsv6MpFpWn4R7byqVyqqDrlKpLOtreHi47/Bciz6vx7Ky7m7ntlKprPq8t46T1Ue7mvp9X+SpfXHcbn3kfe/1cwyDkHUMRRtkfcoZ7m607c52VdJTEXFLm31PSfrdiPjn5vozkn4jIq754hjbB9SYutHU1NTOl1/O9RUJkDQ0NKQ8r1U7tvXGG28sW+/XWvR5PaysW+p8bhePaTXnvXWcrNeuXU159VP74rjd+ui2v13bomQdQ9EGWZ/t0xFRyxyzp17buyxpa8v6ZHPbNSLiaETUIqI2MZH5pWZoMTU1NbDnDg8P91vOmvR5PbQ7j53O7dTU1KrPe+vzsvoY5Gvbuj3vuN36yFtfP8cwCFnHULRC6stzea/u0zLvlfS0JEu6Q9K/5emTaZneMOfe28KcO3PuN0Jdi27IOXdJj0t6RdJ3Jc1Luk/S/ZLub+63pCOSvirpy8ox3x6E+6rU6/WlubvFOe5SqRSlUilsR6VSienp6ahUKkvrnd4809PTYXvpjTY+Pr4UCFljdOuzn7n3jRs3djyWUqkU4+PjS203bdoUmzZtahu609PTy85XVt1Zbev1epRKpaUxSqVSx3PVaZxe2vYqq/Y842adq9b9ed8P11svr3cRBlVf3nDPNee+Fmq1WvB97gDQm+s55w4AuMEQ7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACcoV7rZ32z5ve872oTb7p2w/a/sF2y/afs/gSwUA5JUZ7raHJR2RdI+kHZL2296xotlvSToWEbdK2ifpTwZdKAAgvzxX7rdLmouICxFxRdITkvauaBOS3t58fLOkrw+uRABAr/KE+xZJl1rW55vbWv22pHttz0s6IelD7TqyfcD2rO3ZhYWFVZQLAMhjUDdU90v6VERMSnqPpE/bvqbviDgaEbWIqE1MTAxoaADASnnC/bKkrS3rk81tre6TdEySIuJfJW2UVB5EgQCA3uUJ91OSttveZntUjRumx1e0uShplyTZ/mE1wp15FwAoSGa4R8RVSQ9IOinpnBr/K+aM7Udt72k2e0jSB2x/SdLjkt4fEbFWRQMAuhvJ0ygiTqhxo7R12yMtj89KumuwpQEAVotPqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJChXuNvebfu87Tnbhzq0+TnbZ22fsf0Xgy0TANCLkawGtoclHZH005LmJZ2yfTwizra02S7pYUl3RcQ3bX/vWhUMAMiW58r9dklzEXEhIq5IekLS3hVtPiDpSER8U5Ii4tXBlgkA6EWecN8i6VLL+nxzW6t3Snqn7X+x/Zzt3YMqEADQu8xpmR762S7p3ZImJf2T7R+NiG+1NrJ9QNIBSZqamhrQ0ACAlfJcuV+WtLVlfbK5rdW8pOMR8d2I+Jqkr6gR9stExNGIqEVEbWJiYrU1AwAy5An3U5K2295me1TSPknHV7T5GzWu2mW7rMY0zYUB1gkA6EFmuEfEVUkPSDop6ZykYxFxxvajtvc0m52U9Jrts5KelfRrEfHaWhUNAOjOEVHIwLVaLWZnZwsZGwDWK9unI6KW1Y5PqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AE5Qp327ttn7c9Z/tQl3Y/Yzts1wZXIgCgV5nhbntY0hFJ90jaIWm/7R1t2m2W9KCk5wddJACgN3mu3G+XNBcRFyLiiqQnJO1t0+53JP2epO8MsD4AwCrkCfctki61rM83ty2xfZukrRHxt906sn3A9qzt2YWFhZ6LBQDk0/cNVdtDkj4q6aGsthFxNCJqEVGbmJjod2gAQAd5wv2ypK0t65PNbYs2S7pF0j/a/g9Jd0g6zk1VAChOnnA/JWm77W22RyXtk3R8cWdEvB4R5YioRkRV0nOS9kTE7JpUDADIlBnuEXFV0gOSTko6J+lYRJyx/ajtPWtdIACgdyN5GkXECUknVmx7pEPbd/dfFgCgH3xCFQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQbnC3fZu2+dtz9k+1Gb/r9o+a/tF28/Yrgy+VABAXpnhbntY0hFJ90jaIWm/7R0rmr0gqRYR75L0pKTfH3ShAID88ly53y5pLiIuRMQVSU9I2tvaICKejYhvN1efkzQ52DIBAL3IE+5bJF1qWZ9vbuvkPklP91MUAKA/I4PszPa9kmqSfqrD/gOSDkjS1NTUIIcGALTIc+V+WdLWlvXJ5rZlbN8t6cOS9kTE/7XrKCKORkQtImoTExOrqRcAkEOecD8labvtbbZHJe2TdLy1ge1bJX1CjWB/dfBlAgB6kRnuEXFV0gOSTko6J+lYRJyx/ajtPc1mH5E0LukvbX/R9vEO3QEAroNcc+4RcULSiRXbHml5fPeA6wIA9IFPqAJAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AE5Qp327ttn7c9Z/tQm/1vs/2Z5v7nbVcHXeiimZkZVatVDQ0NqVqtamZmJrN9uVyWbdlWuVy+5jlZfbbuL5fLS/0NDQ0t9dttydtuLftpd9wAEhYRXRdJw5K+KukHJI1K+pKkHSvaHJT0p83H+yR9JqvfnTt3Rq/q9XqMjY2FpKVlbGws6vV6x/ajo6PL2kuKDRs2LD0nq892+9frMjo62vFcAVgfJM1GRr5GRK5wv1PSyZb1hyU9vKLNSUl3Nh+PSPqGJHfrdzXhXqlU2oZWpVLpqX3rc7L67NbHelw6nSsA60PecM8zLbNF0qWW9fnmtrZtIuKqpNcllVZ2ZPuA7VnbswsLCzmGXu7ixYsD2d66L+u53fpYj1I7HgDtXdcbqhFxNCJqEVGbmJjo+flTU1MD2d66L+u53fpYj1I7HgDt5Qn3y5K2tqxPNre1bWN7RNLNkl4bRIGtDh8+rLGxsWXbxsbGdPjw4Y7tR0dHr9m+YcOGpedk9dlu/3o1Ojra8VwBSEzWvI0ac+gXJG3TmzdUf2RFmw9q+Q3VY1n9rmbOPaJxg7NSqYTtqFQqmTcI6/V6lEqlpTnnUql0zXOy+mzdXyqVlvqznWueO2+7teyn3XEDWH+Uc87djbbd2X6PpI+p8T9nHouIw7YfbQ5y3PZGSZ+WdKuk/5K0LyIudOuzVqvF7Oxs5tgAgDfZPh0Rtax2I3k6i4gTkk6s2PZIy+PvSPrZXosEAKwNPqEKAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCcn2IaU0GthckvVzI4O2V1fg2yxRwLDcmjuXGtN6OpRIRmV/OVVi432hsz+b51Nd6wLHcmDiWG1NKx9KKaRkASBDhDgAJItzfdLToAgaIY7kxcSw3ppSOZQlz7gCQIK7cASBBhHsL2x+x/e+2X7T917bfUXRNvbC92/Z523O2DxVdz2rZ3mr7WdtnbZ+x/WDRNfXL9rDtF2w/VXQt/bD9DttPNv+dnLN9Z9E1rZbtX2m+v16y/XjzdymSQbgv93lJt0TEuyR9RdLDBdeTm+1hSUck3SNph6T9tncUW9WqXZX0UETskHSHpA+u42NZ9KCkc0UXMQB/KOnvIuKHJP2Y1ukx2d4i6Zck1SLiFjV+iGhfsVUNFuHeIiI+FxFXm6vPqfF7sevF7ZLmIuJCRFyR9ISkvQXXtCoR8UpEfKH5+H/UCJAtxVa1erYnJb1X0ieLrqUftm+W9JOS/kySIuJKRHyr2Kr6MiLppubvPo9J+nrB9QwU4d7ZL0p6uugierBF0qWW9Xmt40BcZLuqxs83Pl9sJX35mKRfl/RG0YX0aZukBUl/3pxi+qTtTUUXtRoRcVnSH0i6KOkVSa9HxOeKrWqw3nLhbvvvm3NsK5e9LW0+rMbUwExxlcL2uKS/kvTLEfHfRdezGrbfJ+nViDhddC0DMCLpNkkfj4hbJf2vpHV5b8f296jxl+02Sd8vaZPte4utarBy/YZqSiLi7m77bb9f0vsk7Yr19f9EL0va2rI+2dy2LtneoEawz0TEZ4uupw93SdrT/JH5jZLebrseEesxSOYlzUfE4l9RT2qdhrukuyV9LSIWJMn2ZyX9uKR6oVUN0Fvuyr0b27vV+PN5T0R8u+h6enRK0nbb22yPqnFz6HjBNa2Kbasxr3suIj5adD39iIiHI2IyIqpqvCb/sE6DXRHxn5Iu2f7B5qZdks4WWFI/Lkq6w/ZY8/22S+v05nAnb7kr9wx/LOltkj7feL31XETcX2xJ+UTEVdsPSDqpxp3/xyLiTMFlrdZdkn5e0pdtf7G57Tcj4kSBNaHhQ5JmmhcQFyT9QsH1rEpEPG/7SUlfUGMK9gUl9klVPqEKAAliWgYAEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQoP8HD0dAzrFvbm4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y, color='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For the same example above, if we fit a logistic regression model\n",
    "> $ \\displaystyle \\text{logit}(y) = \\log \\left( \\frac{y}{1-y} \\right) = \\beta_0 + \\beta_1 x $\n",
    ">\n",
    "> $ \\displaystyle \\widehat{y} = \\frac{\\exp \\left( \\widehat{\\beta}_0 + \\widehat{\\beta}_1 x \\right)}{1 + \\exp \\left( \\widehat{\\beta}_0 + \\widehat{\\beta}_1 x \\right) } $\n",
    "- The above $\\text{logit}(y)$ is called a link function, there exists other link functions such as probit function and complementary log-log function:\n",
    "> $ \\displaystyle \\text{probit}(y) = \\Phi^{-1}(y) = \\beta_0 + \\beta_1 x $\n",
    "> \n",
    "> $ \\displaystyle \\text{log-log}(y) = \\log(-\\log(1-y)) = \\beta_0 + \\beta_1 x $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=100000.0, class_weight=None, dual=False,\n",
       "          fit_intercept=True, intercept_scaling=1, max_iter=100,\n",
       "          multi_class='warn', n_jobs=None, penalty='l2', random_state=None,\n",
       "          solver='liblinear', tol=0.0001, verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit a logistic regression model\n",
    "clf = LogisticRegression(C=1e5, solver='liblinear') #used for just two classes\n",
    "clf.fit(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Sigmoid(x):\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4FNX6wPHvScEQQtEEC0QTUEBNQg0K8hPxKoIFQQVp0kVJRBAREfGCckFBriKiIChFSQQBUbkWpCqK9BiqNCFAAGlKpEVSzu+PWdLYTTbJzs7u5v08zzyZnZ05553s5s3s2TPnKK01QgghfIuf1QEIIYRwPUnuQgjhgyS5CyGED5LkLoQQPkiSuxBC+CBJ7kII4YMkuQshhA+S5C6EED5IkrsQQvigAKsqDgsL05GRkVZVL0SZdvEi7NwJERFQubLV0Yji2LRp00mtddWi9rMsuUdGRrJx40arqheizPr7b2jWDMaPh4EDrY5GFJdS6oAz+0mzjBBlSFYWdOliJPcBA6yORphJkrsQZciQIZCeDpMmgVJWRyPMZFmzjBDCvT78EL75BtauhcBAq6MRZpPkLkQZsHIlvPIK/PwzXHml1dEId5BmGSF83O7d0KkTzJ0LtWpZHY1wF0nuQviwv/6CNm1g9Gi4+26roxHuJMldCB+VkQHt28ODD0LfvlZHI9xNkrsQPkhrePZZKF/e6M8uyh75QlUIH/Tuu7B6tbH4+1sdjbCCJHchfMx338G4cfDLL1CpktXRCKtIchfCh2zbBj16wJdfggzdVLZJm7sQPuL4caNnzIQJcMcdVkcjrCbJXQgf8M8/8Oij0LWrsQghyV0IL6e10dXxuutg1CiroxGeQtrchfByY8fCjh2wahX4yeWasJHkLoQXW7gQJk82BgMLDrY6GuFJJLkL4aWSkuDpp2HxYqhe3epohKeRD3FCeKEjR6BtW/jgA2jUyOpohCeS5C6Elzl/3kjs/frBY49ZHY3wVJLchfAi2dnQsyfcfDO8/LLV0QhPJm3uQniRV1+Fw4dh+XKZJk8UTpK7EF7i009h9mxYtw6CgqyORng6Se5CeIE1a+C554wr9quvtjoa4Q2kzV0ID3fggPHF6cyZEBNjdTTCW0hyF8KDnTljDAY2ZIgxo5IQzpLkLoSHysqCLl2gSROjSUaI4pDkLoSHGjoUzp2D99+XnjGi+OQLVSE80PTpsGiRMWZMYKDV0QhvJMldCA/zww/GDUo//QRXXWV1NMJbSbOMEB5k717o1Mno0167ttXRCG8myV0ID/HXX/DQQ/Daa3DPPVZHI7xdkcldKTVDKXVcKbXNwfNKKfWuUmqvUmqLUqqh68MUwrdlZMDjj0Pr1sYwvkKUljNX7rOA1oU8fz9Qy7Y8BUwpfVjC0yQmJhIZGYlSioCAgHw/IyMjiY+PJzIyEj8/PyIjI0lMTMx3fHx8PH5+fiilUEpRsWJF4uPjCQsLy9nm7++fU5694y/VV5zlUp1hYWGEhYXh5+dHWFgYFStWdHiMv78/8fHxOefs6Jzy/l7y7pOYmJjvvMLCwkhMTMx3DgEBAcTHxwPGNHkDB8KJE4f54oual/1u7dXrzGtVMG5HcTlz3KXteX+PJYmtpJx5LXy5/pJQWuuid1IqEvhaax1t57mpwA9a6zm2x7uAFlrro4WVGRsbqzdu3FiSmIWbJSYm8tRTT3H+/HmnjwkODmbatGl07dqV+Ph4pkwp3v/8vMcP6tuXDR99RCQQCYQDlYCKQIhtCQD87Sx+tp8l7UmY9zilFFWqVCG4fHkAzl+4wOnTp7H3N+RsfcHBwcwOGMSUs135MrsJIfx9eVlKUbly5Zx6C3P+wgXS0tLyxaSUonz58g5fvypVqgA4PO7ChQt2z7G4sZWUo3Myu17T61+1Cm66qdiHKaU2aa1ji9xRa13kgvE3tc3Bc18D/5fn8XIgtqgyGzVqpIV3iIiI0ECxl4iICK211v7+/sU+9jrQb1aponWjRjrDuLj1yeU7WulrOaL3EWl5LLK4efnttxL9PQIbtS46b7u1K6RS6imMphtuuOEGd1YtSuHgwYOlOi4rK8vpY64F3gCeAAJOn4ZNm3y2v+4ObqE7n/AFj1CDFKvDET7GFX83h4Hr8zwOt227jNZ6GjANjGYZF9Qt3OCGG27gwIEDJToOwN/f36kE3wpIBEILbM8GtgO/AQeAg8CfwFngDHAOyACy7CzZeX6WVN43anj16qxbtw6A22+/ndTDdt/qlx1XUBZhnOJ/VFJDaJY6n9tuu63QssKrV2f9+vVFxuqoHH8/P7Ky7f8Wwm0TsBb3uOLGVlKOzsnsek2vv2rVUkTlBGcu7ym8WeZB4DuMZsYmwHpnypRmGe+RkJCgg4ODNTjfrBIcHKwTEhK01lrHxcUVuX8HuKz55VidOlrPnq2f7927WHW7agkICHB4To5+L4GBgZcdB2ilVJ7H5TSs0jBGx8XFFfk7LlhvcV+r4OBgHRcXpwMDAy8ru1y5cjohIaHQ4wp77YsTW0k5is3sej2l/oJwslnGmcQ+BziKcXGUCvQB+gH9bM8r4H3gd2ArTrS3a0nuXichISGn7f1SG/qlnxERETouLk5HRERopZSOiIi47I0fFxeXL8GFhITouLg4HRoaqpuBvpgnqR/x99fLXnrpsuNL0nZ/qc7Q0FAdGhqqlVI6NDRUh4SEODzGz89Px8XF5Zyzo3PK+3vJu09CQoIODQ3NKS80NFQnJCTouLg47efnr2GGhoW6X7/4In/Hjup15rUqGLejuJw57tL2vL/HksRWUs68Fr5cf17OJnenesuYQXrLCMC4cycqCo7aOlfdcosxI8V111kbl0nefBPmzjWGFqhQwepohDdytreMr35XJbzFSy/lJvawMFi82GcT+5dfwrvvGoOBSWIXZpPkLqyzdStMm5b7eNo08NFeVMnJ0LcvfPsthIdbHY0oC2RsGWGdMWNy1x94ANq1sy4WEx09Cg8/DJMnQ+PGVkcjygpJ7sIa+/fDvHm5j0eN8skZKS5cMP5n9e0LHTpYHY0oSyS5C2vMmGH0jQG47z5o1MjaeEygNfTqZdxh/sorVkcjyhppcxful5UFM2fmPvbRYRBfew0OHICVK33yQ4nwcJLchfutWgWX7virWtUYxNzHzJ0Ls2bBunUQFGR1NKIskuQu3O+rr3LXO3SAcuWsi8UE69bBgAGwbBlcc43V0YiyStrchXtpnT+5t21rXSwmOHgQHn3UmOC6bl2roxFlmSR34V47dkBKirFeqRK0aGFlNC519qzR5fH556FNG6ujEWWdJHfhXitX5q63bOkzTTJZWdC1K8TGGsldCKtJm7twrx9+yF2/+27LwnC1YcMgLQ3mz5eeMcIzSHIX7pOdDT/+mPvYR5pkZs6EL74wxozxkQ8iwgdIchfus2MHnDxprIeFwa23WhuPC/z4ozH22Y8/QmjBWUaEsJC0uQv3yTtrzf/9n9e3X/z+O3TsCImJcPPNVkcjRH6S3IX7bNiQu+7lI2idPm3cezVyJNx7r9XRCHE5Se7CfXwkuWdmGlfsLVtCXJzV0QhhnyR34R7//ANbtuQ+ji1yIhmPNWiQ0aL09ttWRyKEY/KFqnCPrVshI8NYv/FGuPJKa+MpocmTjVkA16yBAPnrER5M3p7CPfJetTdoYF0cpbB0qTHs/OrVULmy1dEIUThJ7sI9tm3LXY+JsS6OEtq5E554AhYsMD54COHppM1duEfe5B4dbV0cJXDqlNEzZtw4uPNOq6MRwjmS3IV7eGlyv3jRGOXxscegZ0+roxHCeZLchfn+/NOYJRqMmSu8pF1Da6Or45VXwhtvWB2NEMUjbe7CfNu3567feiv4+1sXSzG89RZs2gQ//wx+chkkvIwkd2G+Xbty173kPv1Fi2DCBGMwsJAQq6MRovgkuQvz7dmTu167tnVxOGnzZnjySfj6a7j+equjEaJk5MOmMN/u3bnrHp7c//jDmE1p0iS47TaroxGi5CS5C/PlTe61alkXRxEuXIB27aB3b2PsGCG8mSR3Ya6sLNi7N/exhyZ3raFPH6hRA0aMsDoaIUpP2tyFuQ4dMjqLA1xzjcfet/+f/8C+fcYUr14+zLwQgCR3Yba8X6Z66FX7vHkwfTqsWwfly1sdjRCuIcldmGvfvtx1D7x5af16eOYZWLYMrr3W6miEcB1pcxfm2r8/d71GDevisOPQIWNogenToV49q6MRwrUkuQtzeWhyP3vW6PI4cKDxUwhfI8ldmMsDk3t2NnTrZgwr/8ILVkcjhDmkzV2YywOT+8svG2OZffaZ9IwRvkuSuzDP2bNw8qSxXq4cVKtmbTzAxx/D/PlGz5hy5ayORgjzSHIX5sl71R4ZafnQij//DEOGwA8/QFiYpaEIYTppcxfmSUnJXY+MtCoKwOiR2aEDzJ5tjDoshK+T5C7Mk5qau37DDZaFkZYGbdrA8OHQqpVlYQjhVpLchXnyJvfwcEtCyMyETp2gRQvo39+SEISwhCR3YR4PSO6DBxtjl02caEn1QlhGvlAV5rE4uX/wASxZAmvWQIC800UZI295YR4Lk/uyZfDqq0YPmSpV3Fq1EB5Bkrswh9aWJfddu6BrV2O0x5tuclu1QngUaXMX5jh9Gs6fN9ZDQqBSJbdUe+oUPPQQvPEG3HWXW6oUwiNJchfmKHjV7ob7/C9ehPbtc6fKE6Isk+QuzJE3uVevbnp1WhvjsleqBGPHml6dEB5P2tyFOdzc3j5hgjHxxurV4O9venVCeDxJ7sIchw/nrpuc3L/+Gt56y+jyGBJialVCeA1J7sIcbrpy37LFaF9ftMjSEQ6E8DjS5i7M4YbkfuyYMYvSxInQpIkpVQjhtSS5C3OYnNzT0+GRR6BHD+jc2eXFC+H1JLkLc5iY3LWGPn3g+uth5EiXFi2Ez5A2d+F6Z84Y4+wCXHEFhIa6tPgxY2DPHvjxR8vn/xDCY0lyF65XsKeMC29gmj8fpk0zpskrX95lxQrhcyS5C9cz6QamjRshPt4Y6fG661xWrBA+ST7UCtczob398GFjWIFp06BBA5cUKYRPk+QuXM/Fyf3cOaPLY//+Rg8ZIUTRJLkL13Nhcs/Ohu7dIToahg4tZVxClCHS5i5cz4VDD/z733D8OHz6qVsGlhTCZ0hyF67noiv32bNhzhyjZ8wVV7ggLiHKEKeaZZRSrZVSu5RSe5VSL9l5vqdS6oRSKtm2POn6UIXXcEFyX73amNz6f/+DqlVdFJcQZUiRV+5KKX/gfaAlkApsUEot0lrvKLDrZ1rr/ibEKLxJejqcPGmsBwTA1VcXu4iUFGPSjU8+gago14YnRFnhzJX7bcBerfU+rfVFYC7Q1tywhNfK295erVqxB1f/+29jmrxhw6B1axfHJkQZ4kybe3XgUJ7HqcDtdvZ7TCnVHNgNDNJaH7KzT45Tp04xa9YsZ+MU3uL0aejZ01ivVAmK8RprDdu2Qdu2xT5UCFGAq7pC/g+I1FrXBZYCH9vbSSn1lFJqo1JqY0ZGhouqFh7ln39y14v5Lei+fUbXx5tucnFMQpRBSmtd+A5KNQVe1Vq3sj0eBqC1fsPB/v7An1rryoWVGxsbqzdu3FiioIUHGzcOXrJ95z5oELz9tlOHTZtmzKa0di1ceaWJ8Qnh5ZRSm7TWsUXt58yV+wagllKqhlKqHNAJWFSgsrwjfTwM/FacYIUPKUFPmRUrYMQIY7o8SexCuEaRbe5a60ylVH/ge8AfmKG13q6UGgVs1FovAgYopR4GMoE/gZ4mxiw8WTGT++7dxmQbn30GtWqZGJcQZYxTNzFprb8Fvi2wbUSe9WHAMNeGJrxSMe5O/fNPo2fMmDHQooW5YQlR1sjYMsK1nLxyz8iADh2gTRt4Um55E8LlJLkL18nIgD/+MNaVcjjoutbGCI/BwfDmm26MT4gyRMaWEa5z9KiRuQGuvRYCA+3uNnEirFljDDFQzHuchBBOkuQuXMeJGZi+/da4Wl+zBipWdFNcQpRBktyF6xTR3r5tm3Hz6ldfQUSE+8ISoiySNnfhOoUk9+PHjS9PJ0yApk3dHJcQZZAkd+E6DpL7P//Ao49C167GIoQwnyR34Tp2krvW0Lev0XFm1CiL4hKiDJI2d+E6dpL72LGwYwesWgV+cikhhNtIcheuUyC5L1wIkycbg4EFB1sXlhBlkSR34RpZWUY/d5uk4+E8/TQsXuywV6QQwkTyQVm4xvHjkJkJwJEro2j7+BVMnQqNGlkclxBllCR34Rq2JpnzlOfh9M+IizN6yAghrCHJXbhGairZKHrwMbdedYxhMkaoEJaS5C5c4+BBRvIaR6jGh/cvRCmrAxKibJPkLlwicXEoCTzBFzzCFTc6NwOTEMI80ltGlNqaNTBoRRtW0IyrOSEDxwjhAeTKXZTKgQPw2GMw6/p/E812Y2NkpKUxCSEkuYtSOHPGmCbvxRfhgb8Sc5+QK3chLCfJXZRIVpYxsfUdd8DA3meMCVEBypUzJuoQQlhKkrsokRdfhAsX4L33QB08kPvEDTfIIDJCeAD5QlUU20cfwf/+Z4wZExiI0fB+ibS3C+ERJLmLYlm5EoYPh59+gquusm1MScndQZK7EB5BPj8Lp+3ZA506wZw5ULt2nifyXrnLl6lCeARJ7sIpf/1l9IwZNQr+9a8CT8qVuxAeR5K7KFJGBjz+ONx/Pzz9tJ0d5MpdCI8jyV0USmsYOND44vSttxzsJFfuQngc+UJVFOq994wp8n75Bfz97exw/rwxljtAQABUq+bW+IQQ9klyFw4tXgyvv24k9kqVHOx08GDu+vXXO/gPIIRwN0nuwq4dO6B7d/jiC6hRo5Ad8zbJSHu7EB5D2tzFZU6cgDZtjDb2Zs2K2Fna24XwSJLcRT7//GNMj9exI3Tr5sQBe/fmrtesaVpcQojikeQucmhtdHW8+moYPdrJg3bvzl2vU8eUuIQQxSdt7iLHm2/C1q1G7xinx/7Km9zz3bYqhLCSJHcBwJdfwqRJsG4dVKjg5EGZmfD777mPb7rJlNiEEMUnyV3w66/Qty989x1Ur16MA1NSjAQPxoEhIWaEJ4QoAWlzL+OOHoW2bWHKFIiNLebBeZtkatVyaVxCiNKR5F6GXbhgJPannoL27UtQwK5dueuS3IXwKJLcy6jsbOjZ0/gOdPjwEhayfXvuelSUK8ISQriItLmXUa+9BocOwYoVoFQJC9m6NXc9OtolcQkhXEOSexk0Zw58/LHRMyYoqISFZGfnv3KX5C6ER5HkXsasXWsM4bt8OVxzTSkKSkmBc+eM9apVS1mYEMLVpM29DDl40BhaYMYMiIkpZWHbtuWuy1W7EB5HknsZceaMMRjYCy8Y0+WV2ubNueuS3IXwOJLcy4CsLHjiCbjtNhg0yEWFbtiQu96woYsKFUK4irS5lwHDhkFaGsyfX4qeMXlpnT+533abCwoVQriSJHcfN3OmMeHG2rVQrpyLCj18GP74w1gPCZHRIIXwQJLcfdiPP8JLLxk/Q0NdWHDeq/ZGjWRqPSE8kLS5+6jffzcm3EhMhJtvdnHha9bkrjdu7OLChRCuIMndB50+bfSIGTkS7r3XhApWrsxdb97chAqEEKUlyd3HZGYaV+wtW0JcnAkVpKVBUpKx7ucHd95pQiVCiNKS5O5jnnvO6BHz9tsmVfDTT8bQAwANGkCVKiZVJIQoDflC1Ye8/77RYvLLLxBg1iv7/fe56y1amFSJEKK0JLn7iCVLjEmtV6+GypVNqkRrYz6+Sx580KSKhBClJcndB/z2m3EH6uefQ82aJlaUlASpqcb6lVdKe7sQHkza3L3cyZPGmDHjx7sh186bl7vepo2JbT9CiNKS5O7FLl6Exx4zpsjr0cPkyjIyjEHgL+nQweQKhRClIcndS2kN/frBVVfB66+7ocKvv4Zjx4z1atWgdWs3VCqEKCn5XO2l/vtf+PVXo2ein9n/orWGceNyH/fsKU0yQng4+Qv1QosWwcSJxmBgISFuqHDxYmNOPjBGH4uPd0OlQojSkOTuZZKToU8f+PZbCA93Q4Xp6ca8fJf07QvVq7uhYiFEaUibuxf54w9o29a4Wclt43UNGQJ79hjrlSvD8OFuqlgIURqS3L3EhQtGYu/TBx5/3E2VvvMOvPde7uOxY+G669xUuRCiNCS5ewGtoXdvuPFG+Pe/3VDhhQvGfHx55+Rr3x6eftoNlQshXMGp5K6Uaq2U2qWU2quUesnO81copT6zPb9OKRXp6kDLsv/8B/bvh+nTXTRNniPHjhl3Q91yi3HVfkmzZsaUTqZWLoRwpSKTu1LKH3gfuB+4FeislLq1wG59gL+01jcBE4BxCAASExOJjIzETykiIyJITEgwRlW8tGRl8ens2dSMiCBAKa7w9ydAKW6KiODTjz9m3qeZTHr3LMdSY6kSrKgdEcGcWbMgPZ05M2dSJyKCYKUIDw0lrEIFgpXKWSrkWUKU4mqluFEp6inFnUrRWSmGKsVUpdisFNnXXgsvvggHDuTE/xUQsno1qmJFlFLFWsLCwkhMTLTsdy9EWaa01oXvoFRT4FWtdSvb42EAWus38uzzvW2fNUqpAOAPoKoupPDY2Fi9cePG4kWblAR3322s5y3azPXSlFFK62nMg3zDMu6lHltcVq4zTgLDgWmlLKdcuXLMmDGDrl27uiAq82RkZJCamkp6errVoQgBQFBQEOHh4QQGBubbrpTapLWOLep4Z7pCVgcO5XmcCtzuaB+tdaZSKg0IxcgRrpOVBX//7dIiPdUhwnmUhUynj9sSexbwCzAdmA+cd0GZFy9eZPjw4R6f3FNTU6lYsSKRkZEoaX4SFtNac+rUKVJTU6lRo0aJynBrP3el1FPAUwA33HCDO6v2PLYEkq01ea/1NXCWCrRhEf2ZyP38jwzb9rz7FFwv6nkwkvVZ4Izt5x/AAduyHVhr2+5qBw8eNKFU10pPT5fELjyGUorQ0FBOnDhR4jKcSe6HgevzPA63bbO3T6qtWaYycKpgQVrradg+6cfGxha//aJBA2OC0Evy/iG6a70Y+0ZGRnLATmKLiIggJSUFgJqRkRzI08YNCpgN/Mpm/ssr/v5kZWVdVoa/g+2eyFv+kUtiF56ktO9HZ3rLbABqKaVqKKXKAZ2ARQX2WQRcGpewPbCisPb2EgsIMG6kubRUqpS7VKyYu4SE5C4VKuQuwcG5S/nyuUtQUO5yxRW5S7lyuUtgoLEEBOQu/v65i59f7qIUKMWY118nODg43ykEBwczZsyYnMdjxowpsM/rwFVAP4KDg3nqqafslmFvuycqV65cvvMVjoW4YCyJI0eO0L59e4fPnz59msmTJzu9f0E9e/akRo0a1K9fn3r16rF8+fJSxetqH3zwAZ988kmpy0lJSUEpxSuvvJKz7eTJkwQGBtK/f/9ileXM6+qK1/4yWusiF+ABYDfwOzDctm0U8LBtPQijmXYvsB6oWVSZjRo10mVBQkKCjoiI0EopHRERoRMSEhzuAz007NEQmm9fR2Xk3R4aGqorVKigMVpiPGIJDQ21e76eaMeOHVaHoCtUqGB6Hfv379dRUVElPr5Hjx56/vz5WmutV6xYoW+66SaXxJWRkeGSclxl//79ukaNGrp+/fo52yZPnqzr1aunn3nmmWKV5czr6mgfe+9LYKN2Jm87s5MZS1lJ7s5atUrrqlW19oAcUyZ5anLfv3+/vvvuu3VMTIz+17/+pQ8cOKC11nrv3r369ttv19HR0Xr48OE5x+ZN3tu2bdONGzfW9erV0zExMXr37t26Y8eOOigoSNerV0+/8MIL+fbPzMzUgwcP1lFRUTomJka/++67l8WTN7lfuHBBly9fPue5jRs36ubNm+uGDRvq++67Tx85ckRrrfX69et1TExMTp2X6ps5c6Zu06aNvvvuu3Xz5s211lq/+eabOjY2VsfExOgRI0ZorbU+e/asfuCBB3TdunV1VFSUnjt3rtZa66FDh+pbbrlFx8TE6MGDB2uttR45cqQeP3681lrrX3/9Vd9+++06JiZGt2vXTv/5559aa63vuusu/eKLL+rGjRvrWrVq6VWrVtn9vUdFRenOnTvrDRs25Bw3ZsyYnOTu6LXZt2+fbtKkyWWvjaPzc/Taa1265C53qHqAffuMIQUSEoz7h4TFbM1qpizF9Oyzz9KjRw+2bNlC165dGTBgAAADBw5k4MCBbN26lXAHI8h98MEHDBw4kOTkZDZu3Eh4eDhjx47lxhtvJDk5mfHjx+fbf9q0aaSkpJCcnJxTX2EWL15Mu3btAKMr6bPPPsuCBQvYtGkTvXv3ZrhtHKJevXoxdepUkpOT8ff3z1dGUlISCxYs4Mcff2TJkiXs2bOH9evXk5yczKZNm1i1ahWLFy+mWrVqbN68mW3bttG6dWtOnTrFF198wfbt29myZUu+5pNLunfvzrhx49iyZQsxMTG89tprOc9lZmayfv163nnnnXzbC+rUqRNz587l0KFD+Pv7U61atZznCntt4uLi2Lp1K9flGa7D0fmZRZK7xdLS4KGH4JVX4L77rI5GeJo1a9bQpUsXALp168bPP/+cs72DbTasS88X1LRpU15//XXGjRvHgQMHKF++fKF1LVu2jKeffpoA21j9V111ld39hgwZQu3atenSpQtDhw4FYNeuXWzbto2WLVtSv359Ro8eTWpqKqdPn+bMmTM0bdrUbqwtW7bMqWfJkiUsWbKEBg0a0LBhQ3bu3MmePXuIiYlh6dKlDB06lJ9++onKlStTuXJlgoKC6NOnDwsXLrzs+6e0tDROnz7NXXfdBUCPHj3yJdJHH30UgEaNGuV0brCndevWLF26lLlz59KxY8d8zzl6bVavXk3nzp1ztl/i6PzMIkP+WigzEzp2hH/9C555xupohK/p0qULt99+O9988w0PPPAAU6dOpaYLZlAfP3487du3Z9KkSfTu3ZtNmzahtSYqKoo1a9bk2/d03t5tdlSoUCFnXWvNsGHDeNrOGEZJSUl8++23vPLKK9xzzz2MGDGC9evXs3z5chYsWMB7773HihUrnD6HK664AjB6nWVmZjrcr1y5cjRq1IjG3EMRAAAOuUlEQVS33nqLHTt2sGhRwb4k9tnr6VLY+ZlBrtwt9Pzzxg2teYdxER7A+DLKnKWY7rjjDubOnQsYQ1ncaZsFvUmTJnz++ecAOc8XtG/fPmrWrMmAAQNo27YtW7ZsoWLFipw5c8bu/i1btmTq1Kk5ye7PP/8sNLb+/fuTnZ3N999/T506dThx4kROcs/IyGD79u1UqVKFihUrss422YujWAFatWrFjBkzOHvWuNvi8OHDHD9+nCNHjhAcHMwTTzzBkCFDSEpK4uzZs6SlpfHAAw8wYcIENm/enK+sypUrc+WVV/LTTz8BMHv27Jyr+OIaPHgw48aNu+yTjKPXplmzZvm2F3V+ZpErd4tMmQJLl8KaNTJjnTCcP38+X/v5888/z6RJk+jVqxfjx4+natWqzJw5E4B33nmHJ554gjFjxtC6dWsqV658WXnz5s1j9uzZBAYGcu211/Lyyy9z1VVX0axZM6Kjo7n//vt5Js9HxieffJLdu3dTt25dAgMD6du3b6Hd/i51FXzzzTdp1aoVCxYsYMCAAaSlpZGZmclzzz1HVFQU06dPp2/fvvj5+XHXXXfZjRXgvvvu47fffstpwgkJCSEhIYG9e/cyZMgQ/Pz8CAwMZMqUKZw5c4a2bduSnp6O1pq33377svI+/vhj+vXrx/nz56lZs2bO7664oqKiiIqKumy7o9dm4sSJdOnShXHjxtG2bdsiz+/qq68uUVxFKXJsGbOUaGwZH7F0KXTrBqtXG8P4Cuv99ttv3OJF32afP3+e8uXLo5Ri7ty5zJkzh6+++srqsOw6e/ZsTj/usWPHcvToUSZOnGhxVN7B3vvSlWPLCBfauROeeALmz5fELkpu06ZN9O/fH601VapUYcaMGVaH5NA333zDG2+8QWZmJhEREcyaNcvqkMoESe5udOoUtGljTGjUvLnV0Qhvduedd17WzuypOnbseFlPE2E++ULVTS5ehMceg0cegV69rI5GCOHrJLm7gdYQH28Mh/PGG0XvL4QQpSXNMm7w9tuwcSP8/LMxxpgQQphNkrvJvv7aSO5r1hiDVAohhDtIs4yJtmyB3r3h88/BS4Y0FxayN+yrq4awLY4WLVpQp04d6tWrR+PGjUlOTnZr/UUZMWIEy5YtK3U5P/zwA0opPvroo5xtycnJKKX473//63Q5KSkpREdHl3ofV5PkbpJjx+Dhh2HiRGjSxOpohLfq168f3bt3N618rTXZ2dmXbU9MTGTz5s3Ex8czZMgQl9RV2G3+xTFq1Cjuvfdel5QVHR3NvHnzch7PmTOHevXquaRsq0lyN0F6utErpkcPsI0fJESJvPrqqzlXkS1atGDo0KHcdttt1K5dO+fW+qysLIYMGULjxo2pW7cuU6dOBYybh+655x4aNmxITExMzk1OKSkp1KlTh+7duxMdHc2hQ4fsV44x+Njhw7kTry1ZsoSmTZvSsGFDOnTokHMr/bfffsvNN99Mo0aNGDBgAA899FBO/N26daNZs2Z069bNYaxHjx6lefPm1K9fn+joaH766SeysrLo2bMn0dHRxMTEMGHCBMCYMGTBggUALF++nAYNGhATE0Pv3r35559/AGMWtJEjR+ac+86dO+2eX0REBOnp6Rw7dgytNYsXL+b+++/PeT45OZkmTZpQt25dHnnkEf766y/AuM+gXr161KtXj/fffz9nf0fnZwVJ7i6mNfTpA9dfDyNHWh2N8DX2hqqdPn06lStXZsOGDWzYsIEPP/yQ/fv3ExQUxBdffEFSUhIrV65k8ODBlybfYc+ePcTHx7N9+3YiIiIc1pd3WN+TJ08yevRoli1bRlJSErGxsbz99tukp6fz9NNP891337Fp06bL5v3csWMHy5YtY86cOQ5j/fTTT2nVqhXJycls3ryZ+vXrk5yczOHDh9m2bRtbt26lV4E+xOnp6fTs2ZPPPvuMrVu3kpmZyZQpU3KeDwsLIykpibi4uEKbWdq3b8/8+fP55ZdfaNiwYc6gYuB42OBevXoxadKky+41cHR+VpAvVF1szBjYswd+/NGYcU94JzOmU3XFSB/2hqpdsmQJW7ZsybmaTUtLY8+ePYSHh/Pyyy+zatUq/Pz8OHz4MMeOHQOMK9YmhbQXdu3alYsXL3L27NmcNve1a9eyY8cOmjVrBsDFixdp2rQpO3fupGbNmtSoUQOAzp07M23atJyyHn744Zzhhh3F2rhxY3r37k1GRgbt2rWjfv361KxZk3379vHss8/y4IMPcl+BMbF37dpFjRo1qF27NmAM6/v+++/z3HPPXfa7WrhwocNzffzxx+nYsSM7d+6kc+fO/PLLLzmxFRw2uEOHDpw+fZrTp0/T3HYnYrdu3fjuu+8KPb9LMbqTJHcXmj8fpk2DdeuMqVmF97JoyKUi2RuqVmvNpEmTaNWqVb59Z82axYkTJ9i0aROBgYFERkaSnp4O5B9q157ExEQaNWrEkCFDePbZZ1m4cCFaa1q2bMmcOXPy7VvUF64Fh/W1FyvAqlWr+Oabb+jZsyfPP/883bt3Z/PmzXz//fd88MEHzJs3r1jDLDg7rO+1115LYGAgS5cuZeLEiTnJvSQcnV9hY8abRa4tXWTjRuNGpa++gjyTrwhhulatWjFlyhQyMjIA2L17N+fOnSMtLY2rr76awMBAVq5cyYEDB4pVrlKK//znP6xdu5adO3fSpEkTVq9ezd69ewE4d+4cu3fvpk6dOuzbty8ngX322WfFjvXAgQNcc8019O3blyeffJKkpCROnjxJdnY2jz32GKNHjyYpKSlfWXXq1CElJSUnntIM6ztq1CjGjRuXb6YoR8MGV6lShSpVquRMzlFwWF9752cFuXJ3gdRUaNcOPvwQGjSwOhrhrewN+euMJ598kpSUFBo2bIjWmqpVq/Lll1/StWtX2rRpQ0xMDLGxsdx8883Fjql8+fIMHjyY8ePHM336dGbNmkXnzp1zvrgcPXo0tWvXZvLkybRu3ZoKFSrQuHHjYsf6ww8/MH78eAIDAwkJCeGTTz7h8OHD9OrVK6c3zxsFbu8OCgpi5syZdOjQgczMTBo3bky/fv2KfY5gjM1uj6Nhg2fOnEnv3r1RSuVrLnJ0flaQIX9L6dw5uPNOY0Yl24xjwgt525C/nubSsL5aa5555hlq1arFoEGDrA7L65VmyF9plimF7Gzo3h3q1oUXX7Q6GiGs8+GHH1K/fn2ioqJIS0tz21RywjFplimFV16BEyfg00/N6V0hhLcYNGiQXKl7GEnuJfTJJ/DZZ0bPmDzdYoUQwiNIci+Bn3+GF16AH36AsDCroxGuorW2O2u9EFYo7feh0uZeTPv3Q4cOxpX7rbdaHY1wlaCgIE6dOlXqPyghXEFrzalTpwgKCipxGXLlXgx//21Mk/fyy9C6tdXRCFcKDw8nNTX1slvnhbBKUFBQvq6xxSXJ3UmZmdCpkzH3af/+VkcjXC0wMDDn9nkhfIE0yzjphRcgI8MYwleaZYUQnk6u3J0wdSosXgxr10JgoNXRCCFE0SS5F2HFChgxAlavhipVrI5GCCGcY9nwA0qpE0DxRjIyVxhw0uogCuHp8YHnx+jp8YHnx+jp8YHvxxihta5a1E6WJXdPo5Ta6Mx4DVbx9PjA82P09PjA82P09PhAYrxEvlAVQggfJMldCCF8kCT3XNOK3sVSnh4feH6Mnh4feH6Mnh4fSIyAtLkLIYRPkit3IYTwQZLcC1BKDVZKaaWUx433qJQar5TaqZTaopT6QinlET3vlVKtlVK7lFJ7lVIvWR1PQUqp65VSK5VSO5RS25VSA62OyR6llL9S6lel1NdWx2KPUqqKUmqB7T34m1KqqdUxFaSUGmR7jbcppeYopUo+8pZr4pmhlDqulNqWZ9tVSqmlSqk9tp9XmlG3JPc8lFLXA/cBB62OxYGlQLTWui6wGxhmcTwopfyB94H7gVuBzkopTxsvMxMYrLW+FWgCPOOBMQIMBH6zOohCTAQWa61vBurhYbEqpaoDA4BYrXU04A90sjYqZgEFhxl8CViuta4FLLc9djlJ7vlNAF4EPPKLCK31Eq11pu3hWqDkQ8a5zm3AXq31Pq31RWAu0NbimPLRWh/VWifZ1s9gJKXq1kaVn1IqHHgQ+MjqWOxRSlUGmgPTAbTWF7XWp62Nyq4AoLxSKgAIBo5YGYzWehXwZ4HNbYGPbesfA+3MqFuSu41Sqi1wWGu92epYnNQb+M7qIDCS5KE8j1PxsMSZl1IqEmgArLM2ksu8g3FhkW11IA7UAE4AM21NRx8ppSpYHVReWuvDwH8xPnkfBdK01kusjcqua7TWR23rfwDXmFFJmUruSqlltra4gktb4GVghIfHeGmf4RhNDYnWRep9lFIhwOfAc1rrv62O5xKl1EPAca31JqtjKUQA0BCYorVuAJzDpOaEkrK1XbfF+EdUDaiglHrC2qgKp43uiqa0FJSpgcO01vfa266UisF4Q2y2TbMWDiQppW7TWv/hxhAdxniJUqon8BBwj/aMfqyHgevzPA63bfMoSqlAjMSeqLVeaHU8BTQDHlZKPQAEAZWUUglaa09KTKlAqtb60ieeBXhYcgfuBfZrrU8AKKUWAncACZZGdbljSqnrtNZHlVLXAcfNqKRMXbk7orXeqrW+WmsdqbWOxHgjN3R3Yi+KUqo1xkf3h7XW562Ox2YDUEspVUMpVQ7jC6xFFseUjzL+Y08HftNav211PAVprYdprcNt771OwAoPS+zY/hYOKaXq2DbdA+ywMCR7DgJNlFLBttf8HjzsS1+bRUAP23oP4CszKilTV+4+4D3gCmCp7RPGWq11PysD0lpnKqX6A99j9E6YobXebmVMdjQDugFblVLJtm0va62/tTAmb/QskGj7J74P6GVxPPlordcppRYASRjNlr9i8d2qSqk5QAsgTCmVCowExgLzlFJ9MEbGfdyUuj3jk70QQghXkmYZIYTwQZLchRDCB0lyF0IIHyTJXQghfJAkdyGE8EGS3IUQwgdJchdCCB8kyV0IIXzQ/wNBeqQZau1LzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Predict\n",
    "prob = Sigmoid(clf.intercept_ + clf.coef_ * X_test).ravel()\n",
    "\n",
    "# Plot the result\n",
    "plt.plot(X_test, prob, color='red', linewidth=3)\n",
    "plt.plot(X_test, b0 + b1 * X_test, color='blue', linewidth=1);\n",
    "plt.scatter(x, y, color='black');\n",
    "plt.axhline(0.5, color='0.5');\n",
    "plt.ylim(-0.25, 1.25);\n",
    "plt.yticks([0, 0.5, 1]);\n",
    "plt.legend(('Logistic Regression Model', 'Linear Regression Model'), loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Note:** For `LogisticRegression` in scikit-learn, it minimizes a cost function with some regularization\n",
    "> $ \\displaystyle \\min_{\\omega, c} \\frac{1}{2} \\omega^\\top \\omega + C \\sum_{i=1}^{n} \\log \\left[ \\exp \\left( -y_i \\left( X_i^\\top \\omega + c \\right) \\right) + 1 \\right] \\qquad L_2 \\text{-regularization} $ \n",
    ">\n",
    "> $ \\displaystyle \\min_{\\omega, c} \\lvert\\lvert \\omega \\rvert\\rvert_1 + C \\sum_{i=1}^{n} \\log \\left[ \\exp \\left( -y_i \\left( X_i^\\top \\omega + c \\right) \\right) + 1 \\right] \\qquad L_1 \\text{-regularization} $ \n",
    "\n",
    "[References](http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression (Multinomial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- General Form (K class)\n",
    "> $ \\displaystyle \\log \\left[ \\frac{\\mathbb{P}(y=1 \\ | \\ X = x) }{\\mathbb{P}(y=K \\ | \\ X = x) } \\right] = \\beta_{10} + \\beta_{11} x $\n",
    ">\n",
    "> $ \\displaystyle \\log \\left[ \\frac{\\mathbb{P}(y=2 \\ | \\ X = x) }{\\mathbb{P}(y=K \\ | \\ X = x) } \\right] = \\beta_{20} + \\beta_{21} x $\n",
    ">\n",
    "> ...\n",
    ">\n",
    "> $ \\displaystyle \\log \\left[ \\frac{\\mathbb{P}(y=K-1 \\ | \\ X = x) }{\\mathbb{P}(y=K \\ | \\ X = x) } \\right] = \\beta_{(K-1)0} + \\beta_{(K-1)1} x $\n",
    ">\n",
    "\n",
    "- Fit (maximum likelihood)\n",
    "> $ \\displaystyle l(\\beta) = \\sum_{i=1}^{n} \\log \\pi_{g_i}(x_i; \\beta) \\qquad $ where $ \\displaystyle \\qquad \\pi_{k}(x_i; \\beta) = \\mathbb{P}(y=k \\ | \\ X = x_i; \\  \\beta) $\n",
    ">\n",
    "> $ \\displaystyle \\widehat{\\beta} = \\arg \\max_{\\beta} l(\\beta) $\n",
    "\n",
    "- Predict\n",
    "> $ \\displaystyle \\mathbb{P}(y=k \\ | \\ X = x) = \\frac{\\exp (\\widehat{\\beta}_{k0} + \\widehat{\\beta}_{k1} x )}{1 + \\sum_{j=1}^{K-1} \\exp \\left( \\widehat{\\beta}_{j0} + \\widehat{\\beta}_{j1} x \\right) } \\qquad k = 1,2, \\cdots, K-1 $\n",
    ">\n",
    "> $ \\displaystyle \\mathbb{P}(y=K \\ | \\ X = x) = \\frac{1}{1 + \\sum_{j=1}^{K-1} \\exp \\left( \\widehat{\\beta}_{j0} + \\widehat{\\beta}_{j1} x \\right) } $\n",
    "\n",
    "- Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 1, 2]), array([50, 50, 50]))"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iris = datasets.load_iris()\n",
    "X = iris.data[:,:2]             # first two features\n",
    "Y = iris.target\n",
    "np.unique(Y, return_counts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=100000.0, class_weight=None, dual=False,\n",
       "          fit_intercept=True, intercept_scaling=1, max_iter=100,\n",
       "          multi_class='multinomial', n_jobs=None, penalty='l2',\n",
       "          random_state=None, solver='lbfgs', tol=0.0001, verbose=0,\n",
       "          warm_start=False)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit the model\n",
    "model1 = LogisticRegression(C=1e5, solver='lbfgs',  multi_class='multinomial') #used for multi-class\n",
    "model1.fit(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(171, 231)"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate data for decision boundary\n",
    "h = 0.02\n",
    "x_min, x_max = X[:,0].min() - 0.5, X[:,0].max() + 0.5\n",
    "y_min, y_max = X[:,1].min() - 0.5, X[:,1].max() + 0.5\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "\n",
    "# Predict\n",
    "pred = model1.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "pred = pred.reshape(xx.shape)\n",
    "pred.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7MAAAIaCAYAAADyehr9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4VNXaxuFnTXpPCBBCDaFKF+mICEhREVRU8AgqHgsesRw7NmzYPhWxC6LHgthRRBQ7SBEBKVKl10CA9F5mfX/A4RBpmZDJzoTffV1zmVnZWTzDSGbeWe9e21hrBQAAAACAL3E5HQAAAAAAAE9RzAIAAAAAfA7FLAAAAADA51DMAgAAAAB8DsUsAAAAAMDnUMwCAAAAAHwOxSwAAAAAwOdQzAIAAAAAfA7FLAAAAADA51DMAgAAAAB8jr/TATwVEV3N1qhd1+kYAAAAlY5rz3qnIwDASduYkrfPWlvjRMf5XDFbo3ZdjZsy0+kYAAAAlU7Y//V1OgIAnLTBU9duLc1xtBkDAAAAAHyOz63MAgAA4ABWYgGcyliZBQAAAAD4HIpZAAAAAIDPoc0YAADAR9BWDAD/w8osAAAAAMDnUMwCAAAAAHwObcYAAACVFG3FAHBsrMwCAAAAAHwOxSwAAAAAwOdQzAIAAAAAfA7FLAAAAADA57ABFAAAQCXBhk8AUHqszAIAAAAAfA7FLAAAAADA59BmDAAA4BDaigGg7FiZBQAAAAD4HIpZAAAAAIDPoc0YAACgAtFaDADlg5VZAAAAAIDPoZgFAAAAAPgc2owBAAC8iLZiAPAOVmYBAAAAAD6HYhYAAAAA4HNoMwYAAChHtBUDQMVgZRYAAAAA4HMoZgEAAAAAPodiFgAAAADgcyhmAQAAAAA+h2IWAAAAAOBz2M0YAADgJLB7MQA4g5VZAAAAAIDPoZgFAAAAAPgc2owBAAA8QFsxAFQOrMwCAAAAAHwOK7MAAADHwUosAFROrMwCAAAAAHwOxSwAAAAAwOfQZgwAAHAY2ooBwDewMgsAAAAA8DkUswAAAAAAn0ObMQAAOKXRVgwAvomVWQAAAACAz6GYBQAAAAD4HIpZAAAAAIDPoZgFAAAAAPgcilkAAAAAgM9hN2MAAHDKYQdjAPB9rMwCAAAAAHwOxSwAAAAAwOfQZgwAAKo82ooBoOrx+sqsMcbPGLPUGDPjKN+72hiz1xiz7ODtWm/nAQAAAAD4vopYmb1V0hpJkcf4/kfW2tEVkAMAAAAAUEV4tZg1xtSVdL6kcZJu9+afBQAA8F+0FQNA1eftNuMXJN0tyX2cY4YYY1YYYz41xtTzch4AAAAAQBXgtWLWGDNQUrK1dslxDvtKUoK1to2k7yW9c4y5rjfGLDbGLM5MTfFCWgAAAACAL/Fmm3F3SYOMMedJCpYUaYx531o7/L8HWGv3H3b8m5KeOdpE1tqJkiZKUmKLNtZ7kQEAgC+irRgATj1eW5m11o6x1ta11iZIGibpp8MLWUkyxsQfdneQDmwUBQAAAADAcVX4dWaNMY9KWmytnS7pFmPMIElFklIkXV3ReQAAAAAAvqdCillr7S+Sfjn49UOHjY+RNKYiMgAAAAAAqg5v72YMAAAAAEC5q/A2YwAAgJPFhk8AAFZmAQAAAAA+h2IWAAAAAOBzaDMGAACVHm3FAIC/Y2UWAAAAAOBzWJkFAHhV6t7dmjrhSS36+Vu5/PzUpe9ADbt5jCKiY5yOBgAAfBjFLADAawry8/T49UPVsfe5mjBjvoqLivTF5Jf05L+u0OPvfSWXn5/TEVFJ0VYMADgR2owBAF7z+w8zVT2+robdfK8iY2IVUyNOV9/zmFwuo+ULZjsdDwAA+DCKWQCA1+zY9Jeat+9UYswYo2and9LOTX85lAoAAFQFtBkDALymTsMmmv/tFyXGrLX6a/liXXzdbQ6lQmVEWzEAwFOszAIAvKbzOedpz46t+viVZ5SVnqr0lH1699mHVVRYqLbdznY6HgAA8GEUswAArwkMDtEDEz9S8s5tumlAR912QXflZWdpzKsfsPkTAAA4KbQZAwC8qlrNeI1+4mVZayUdOGcWoK0YAHCyKGYBABWCIhYAAJQn2owBAAAAAD6HYhYAAAAA4HMoZgEAAAAAPodiFgAAAADgc9gACgAAVAh2MAYAlCdWZgEAAAAAPodiFgAAAADgc2gzBgAAXkFbMQDAm1iZBQAAAAD4HFZmAQBAuWAlFgBQkViZBQAAAAD4HIpZAAAAAIDPoc0YAACUCW3FAAAnsTILAAAAAPA5FLMAAAAAAJ9DmzEAACgV2ooBAJUJK7MAAAAAAJ9DMQsAAAAA8DkUswAAAAAAn0MxCwBVyO5tm7Vl3SoVFxU5HQUAAMCr2AAKAKqA5J3b9OoDtyp55zaFRUQpNztTI8c8oTN6smEPAAComihmAcDHud1uPXvbSJ11waU674rr5PLz07plizT+jusUP/kz1U5o5HRE+Ch2LwYAVGa0GQOAj1u3dKFcLj+dP+IGufz8JEnN2nXU2RcO0y9ffOhwOgAAAO+gmAUAH5e+f5/i6jaQMabEeFy9BKWn7HMoFQAAgHfRZgwAPq5JmzM0+YkxyspIU3hktCTJWqvff/haHXuf63A6+BLaigEAvoRiFgB8XGyt2up14eV6/LrLNPia0QqPjtHP06YqI3W/zjzvIqfjAQAAeAXFLABUAZffep8antZac2Z8qvycbLXt3kvXPvC0AoNDnI4GAADgFRSzAFAFGGPUtf8gde0/yOko8CG0FQMAfBkbQAEAAAAAfA7FLAAAAADA59BmDADAKYTWYgBAVcHKLAAAAADA57AyCwBAFcZKLACgqmJlFgAAAADgc1iZBVCl5eVka+aUSfpjzg/yDwhUt/6D1GfIcPn58+sPAADAl/FuDkCVVVRYoCdu/Idia9XW8H8/qPy8XH351sva8OdS/evxCU7HAwAAwEmgmAVQZS366Vv5+/vrlqdelTFGknTaGV10+4Vnaetfq9WgaQuHEwIAAKCsOGcWQJW1fsUfOuPsfocKWUkKDApW225na/2KJQ4mAwAAwMliZRZAlRVTI067tmw6YnzX5g1qfxY7vKJqYvdiAMCpgpVZAFVWj4FDtPjnb7Xop29krVVxUZG++WCy0vYlq223s52OBwAAgJPAyiyAKiu6ek3dMX6yJj12j9555iEVFRYqvkGi7nn5PfkHBDgdDwAAACeBYhZAlda0bQc988kP2rN9iwICgxRbq7bTkYByRVsxAOBURTELoMozxqhW/YZOxwAAAEA54pxZAAAAAIDPYWUWAAAfQlsxAAAHsDILAAAAAPA5FLMAAAAAAJ9DmzEAAJUYbcUAAByd11dmjTF+xpilxpgZR/lekDHmI2PMBmPMQmNMgrfzAAB8W1FhoZbM/k6zPvyP1q/4Q9ZapyMBAAAHVMTK7K2S1kiKPMr3/ikp1Vrb2BgzTNLTkoZWQCYAgA/au2u7nrppuCJiYlWvcTN988Ek1WvUTLc8/ZoCAoOcjgcAACqQV1dmjTF1JZ0v6c1jHDJY0jsHv/5UUh9jjPFmJgCA73rz8XvVc9Blevitz/XP+57Us5/9InexW99MOdbLDAAAqKq83Wb8gqS7JbmP8f06krZLkrW2SFK6pFgvZwIA+KDMtFRtWLlU515x7aEx/4AADf7naM3/9ksHkwEAACd4rZg1xgyUlGytXVIOc11vjFlsjFmcmZpSDukAAL7GXVwkl8slP7+SZ8gEBAapuKjIoVQAAMAp3jxntrukQcaY8yQFS4o0xrxvrR1+2DE7JdWTtMMY4y8pStL+v09krZ0oaaIkJbZow04fAHAKioqtofj6iZo783OddcGlkiRrrWZ9+JbOOLufw+nKD7sXAwBQOl4rZq21YySNkSRjzNmS7vxbIStJ0yVdJWmBpEsk/WTZlhIAcAzX3P+knh49QisXzlW9xs21bN5PysvJ1pV3PuJ0NAAAUMEq/DqzxphHJS221k6XNFnSe8aYDZJSJA2r6DwAAN+R0KylnvnkR837Zpr2Je1Qn0uGq1Pvc+UfEOh0tDJjJRYAgLKpkGLWWvuLpF8Ofv3QYeN5ki6tiAwAgKohIjpGAy6/xukYAADAYd7ezRgAAAAAgHJX4W3GAACcymgrBgCgfLAyCwAAAADwORSzAAAAAACfQ5sxAABeNjT92UNfz3AwBwAAVQkrswDgsJysDBUVFTkdAwAAwKewMgsADvn2g8n6YvJLykxPlZ+/vxKatdK9r7yn0PBIp6MBAABUehSzAOCART99o49eeUajHnleHXsNUErybr35+N0ae/WF+r9Pf3I6Hk7S4W3FAADAO2gzBgAHfPr68xo88iZ1Pud8ufz8VD2+jm595g3t3bVDW9aucjoeAABApUcxCwAOyM5IU+PW7UuMhYSFq2btetqwcqlDqQAAAHwHbcYA4IDImOpa+ftctep85qGxzLRU7dm5TS06dHUwGcqCtmIAACoeK7MA4IAr/v2Avp36lmZ9+LYyUvdr46plemr0cNVv3Ey1Exo5HQ8AAKDSY2UWABzQslN33fjIeH0w4XF98MLjCggMVvP2nXXbM687HQ0AAMAnUMwCgEM69z1fnfue73QMAAAAn0SbMQAAAADA51DMAgAAAAB8Dm3GAAB46GR2Lx54ff8S92dMnHWycQAAOCWxMgsAAAAA8DkUswAAAAAAn0ObMQAAJ3AybcUAAMA7KGYBVDpr/lioSY/epdzsLAWHhumf9z+pVp3OdDoWyshaq6W//qhFP30jl5+fuvS7QK06nSljjNPRAAAVbN2+XP28OV15RW6dUTtc3epFyM/l/OuB21ot3pWl37ZnyeWSetSPVJu4UF6rKjmKWQCVyo+ff6D3/m+suvS7QC07dteaP37Ts7eO1LBbx2jAsGucjgcPWWv15uP3aMOff6jPkOEqLi7W20/er469BujyW+9zOt4xsRILAOVv+toUfbEuRec3iVFEoJ++Wpei2VvSNaZHXUcLWmutXvl9t9bvz1P/xtEqdlu9tmi3utWL0JXtajqWCydGMQugUvn45ad1yY13aOCVoyRJPQYOUf3GzTXtjRcoZn3Q+hVLtGrRPD310fcKDgmVdOA5vWtIb511waWqk9jE4YQAgIqQllekD1fu04RzG6pGWIAkqXdilO79Yat+25Gp7vUjHcu2em+uViXnaMK5DRXkf2BLoV4No3TT15vUOzFKdSODHMuG42MDKACVSnZmhnpddHmJsV4X/UNZGekqLi52KBXKatm8n9Wt/+BDhawkhUdGq2PvAVo+/xfHcgEAKtbK5By1rBl6qJCVJH+XUe+GUfojKdvBZNIfSdnq0SDyUCErSRFBfupSL0JLHc6G42NlFkCl4u/vr8zUFIVFRB0ay0xLkZ+/v/z8/BxMhrIICQ1X8s5tR4xnpqWq4WltHEh0dLQVA4B3hfi7lJl/5IfSmfnFCvF3dn0tJMClfTmFR4xn5BcrJIC1v8qMZwdApRIRE6v3xz+mosIDLyrFRUV6//lHFRkT63AylEW3AYO18IcZ2rJ25aGxtX8s1Krf56pj7wEOJgMAVKQ2cWFKzi7Ugu2Zh8Z2ZxXomw1p6tUw6jg/6X1nNYjU3K0Z2pSad2hsZXKOVu7JUZe6EQ4mw4mwMgugUnnkP9N079D+urHv6WrU6nRtWrVM1krjpsx0OhrKILZWbV37wDN64sbLldC8tYqLCrVz03rd/OQrCo+MdjoeAKCCBPgZjelRR0/N3alpa/YrPNBPa/fl6sp2NdSoWrCj2WqGBejGjrX00E/blBgTrGJrtT29QHd2r63wQLrCKjNjrXU6g0cSW7SxvKkFqr5vP3xLK+bPVqtOZ+q84dc5HQcnKS8nW6sWzZPL5aeWnborMMjZNy6Vqa14xsRZTkcAgApT5LZauSdHuUVutY4LrVTFYm6hWyv2ZMvPZdQmLlSBfjSxOmXw1LVLrLUdTnQcK7MAKqUBw65h9+IqJDg0TGf07Od0DACAw/xdRu3iw5yOcVQhAS51pq3Yp/BxAwAAAADA57AyCwA4JVSm1uLDDby+f4n7tB0DAFA6rMwCAAAAAHwOxSwAAAAAwOdQzAIAAAAAfA7FLIBKKSU5SasWzdf+3bucjnKEpK2btHrxAmWlp570XDlZmVq9eIF2bFxXDskAAABOHWwABaBSKSos1NtP3a/ff5ypeo2ba8fGdWp/1jm69oGn5R8Q6Gi2zLRUvXzfaG3fsFZxdRto+8a/dN4V1+qi626VMcbj+b75YLI+e2O86iY20f49uxRbq45uffpVxdSo5YX0AAAAVQvFLIBK5YvJL2n/7l168evfFBIWrrzcHL1070369PXnNezmex3NNumxu1S7YWPdNeE/8g8IUOrePXryX/9QfINEde0/yKO5ViyYrW+nvqUnPpipmnXqy11crM8mjtfL992sByd94qVHcGqprLsXAwCA8kGbMYBK5edpH2j4HQ8pJCxckhQcEqoRdzykn6dNdTRXeso+rV78m4aNvlf+AQGSpJgacbpk1J36qQzZfpo2VRdeM1o169SXJLn8/HTxdbcpaesmJW3bXK7ZAQAAqiKKWQCVSlZGmmLjapcYi61VW1kZabLWOpRKysnMUFhklIJCQkqMV4urVaZzZ7PT01QtLr7EmJ+/v6Krxykr7eTPxQUAAKjqaDMGUKm06NBN87/5Qn0uGX5obN43X6hFh25lOi+1vBxoBS7SX8uXqGnbMw6Nz//mS7Xo0NXj+Vp07Kr5336ptt3OPjS2c9N67d+9Uw2anlYekU85tBUDAHBqoZgFUKkMHX2PnvrXFdqbtEPNTu+k9SuW6KfPpuiuF//jaC4/f38Nv/0hvXDX9Tpv+HWqndBYS2Z/p5UL5+rht6d5PF/fy67SIyMv1sv336yu/S7Q3l07NOPd1zXsljEKDA458QQAAACnOONk215ZJLZoY8dNmel0DABetGfHVn330X+0c9N61U5orH5Dr1Kt+g2djiVJ2rR6uX749H2lJu9W49bt1feyKxUZE1umubIz0/Xjp+9rzZLfFFmtunpddLman96pnBOfOqrKyuyMibOcjgAAgKMGT127xFrb4UTHUcwCAHxSVSle/45iFgBwqittMcsGUAAAAAAAn8M5swAAn1BVV2L/buD1/UvcZ6UWAICjY2UWAAAAAOBzKGYBAAAAAD6HNmMAQKV0qrQVAwCAsmFlFgAAAADgc1iZBf5m1e/zNOPd15W0bZPqJjbVoJE3qWnbE+4M7nN2bdmoLya/pPXLFyuqek31u+xKdRtwYZnmKios1LdTJ2vezGkqKMhX+x59NGjkaEVEx5RzagAAyqaw2OqrdSmaszVDBcVWHeuE65IWsYoI8nM6GoAyopgFDrNk9veaPO5eXX7LfWrcpr3W/rFQz99xnW556lW16NDV6XjlZs/2LXrsuks14PJrdOE/b9bubZs1dcITSknerYFXjvJ4vlcfvFXZGWm66p7HFBIWrh8+eU+PXXeJHnt3hoJCQrzwCFAV0VYMwJvGL9ilnEK3bugQp2B/l75Zn6b7f9ym/+vXQEH+NCsCvohiFjjMp689qxvGPqu23XtJkuLrN1RQSIg+e+N5tejwicPpys/X709U74v/ocHXjJYk1U5opLqNmuqB4eer76VXeVSAbvtrjf5atkjjp89VQGCQJOma+57QM7dcpQWzvtTZFw7zymMAAKC0NqXmad3+XL0+sJEC/Iwk6caOcXrklx2auy1DfRKjHU4IoCz4GAo4qKiwUNs2rFXrrj1LjLfr3lubVi13KJV3bFq1XO269y4xVrNOfUVXr6nd2zd7Ntfq5WrZ6cxDhawkGWPUrntvbVxdtf7eAAC+aUNKntrEhR0qZKUDr1Vn1A7T+pQ8B5MBOBkUs8BBfv7+io6toZ2b/ioxvm39GlWPr+tQKu+oHl9X29avKTGWk5WplOQ9iqkR5/Fc2zeslbW2xPj2DWtVo4r9vQEAfFPN0ABtTc8/YnxLWr5qhgY4kAhAeaCYBQ4yxmjAP67VpEfvVvLObZIObJL09pP369wrrnU4XfkacPk1+nzieK1btkiSlJG6X5Meu1sdzu6nyJhYj+Zq0bGb3MXF+vT155Sfmyt3cbHmfTNNi3+ZpbMuuNQb8QEA8EjruFAVFVtN/XOv8ovcKnZb/bIlXYt2Zql3wyin4wEoI86ZBQ5z/ojrVZifpwdGDJS/f4DcbrcuuGqUel10udPRylXz9p014o6xevXB25Sfm6PCgnx16z9YI+54yOO5XC6X7nn5Xb35+Bj9q197ufz8FFcvQXe+8Jaiq9f0QnoAADzj5zIae3Zdvbpot66atkF+Lik+IlAP9qyr6BDeDgO+yvy9NbCyS2zRxo6bMtPpGKjiCgvylZGaoqhqsfIPCHQ6jte43W6l7duj0IgoBYeEnvR8WRlpKi4sVFRsjXJIh6qO3YvLZsbEWU5HAHxaVkGxiootRSxQiQ2eunaJtfaE18bkXzFwFAGBQYqNi3c6hte5XC5Vq1l+jzM8kt0gAQCVW3gg15UFqgrOmQUAAAAA+BxWZgEAFYK2YgAAUJ5YmQUAAAAA+ByKWQAAAACAz6HNGADgNbQWAwAAb/FaMWuMCZY0R1LQwT/nU2vt2L8dc7Wk/5O08+DQy9baN72VCYD3bP1rtaa9+aIK8nJ1ziUj1P6sc8o8V0Fenj5/8wVtWrVCCc1a6uIb/n1Slw7atWWj/pjzvfwDAtX5nPMUU6NWmeeqzDJS92vh9zOUm5Otdt17qX6T05yOBAAA4DXebDPOl9TbWttWUjtJA4wxXY5y3EfW2nYHbxSygA9677lHNPaqwfL3D1D1+Dp65YFb9Oh1l5Zprp2b1+umAR20fN7PSmjWUqsWzdPo/h21fcPaMs33+aQX9Oi1l2hf0k5t+2u17r70HM3/9osyzVWZ/THnB9158dla/+dSpe1N1tOjR+j95x+Vr11LHAAAoLS8tjJrD7yDyjp4N+DgjXdVQBWzd9d2/fDJe3r03elq0LSFJOmSUXfq7kt668fPpqjPkCs8mm/C3aPUbcCFuvqex2SMkXSgWJ5w9yg9+/kvHs21afUK/fTZFD390XeKiq0hSTpv+HV6+Johat2lpyKiYzyar7LKy83R6w/frrtffFeNW58uSRoy6t966MpBatu9l1p37lFhWWgrBgAAFcWrG0AZY/yMMcskJUv63lq78CiHDTHGrDDGfGqMqefNPADK39fvvaFmp3c6VMhKUmRMNfUberV+njbV4/l2b9uiC/85+lAhK0kX/vMW7d6+VW6326O5Fv7wtXoOHnqokJWkuo2aqXWXHvpjzvceZ6usVv0+Vw2atjxUyEpSWESU+gwZroXff+1gMgAAAO/xajFrrS221raTVFdSJ2NMq78d8pWkBGttG0nfS3rnaPMYY643xiw2xizOTE3xZmQAHrLWyuU68lfJgbGyNWMY4/rbfZVtLmuPmEuSjFSl2m8PPAfmiHHjclWpxwkAAHC4Crk0j7U2TdLPkgb8bXy/tTb/4N03JZ1xjJ+faK3tYK3tEBFTzbthAXjkvBHXa82S30qc05qVnqpZH/1HZw26zOP54uo20PT/vFpibPrbr6pm3QZHLZqPp2OfczV7+kfKOOxDsJ2bN2jFb7+e1AZVlU2rTmdq85o/tWn1ikNjOVmZ+vHT99Wpz7le/bOHpj9b4gYAAFBRvLmbcQ1JhdbaNGNMiKS+kp7+2zHx1tqkg3cHSVrjrTwAvCOuTgP1umiYHrxykLr2u0ChEZGa89UnqpPQWP0uu8rj+W555nU9cs1FWr9iidp06amVv8/Vjo1/6YFJH3k8V+NWp+usCy7VvUP7qmv/QSrIy9PCH2boyrseVmRMrMfzVVbBoWG67qH/01M3XaEOZ/dXeFSMFnw3XWf07Kc2XXs6HQ8AAMArzIla0IwxHST1kFRbUq6klTpw/mvqCX6ujQ60DfvpwArwx9baR40xj0pabK2dbox5UgeK2CJJKZJutNYed8vSxBZt7LgpM0v14ABUnA0rl+mLyS+qIDdXfYYMV+e+55d5rrycLH3y2nPavGalGjRroUtvvEOh4ZFlnm/7hrVaMvt7+QcEqEvfC1Q9vk6Z56rMUvfu0W/ffaW83Gy17dZLiS3alPufwepr5TNj4iynIwAAUK4GT127xFrb4UTHHbOYNcaMlHSzpM2SlujAJk7BkppK6q4DRe2D1tpt5RW6NChmAcA5FLOVD8UsAKCqKW0xe7w241BJ3a21uUf7pjGmnaQmkiq0mAUAAAAA4JjFrLX2leP9oLV2WfnHAQAAAADgxE64AZQxpqEOtBsnHH68tXaQ92IBAAAAAHBspdnN+AtJk3XgmrBu78YBAAAAAODESlPM5llrX/R6EgBApcOGTwAAoLJyleKYCcaYscaYrsaY9v+9eT0Z4KCiwkKlp+xTcVGR01GOkJWRrp2b18vtrnyNEmn7krVn59Zymau4qKhcn4OsjDTlZmeVy1zlye12a+fm9cpKP+7VzuBDcvIKtT89Vye69J0TsguKlV1QXC5zFRZbpeUVqdhd+R4nAODUUJqV2daSRkjqrf+1GduD94EqxVqrr/7zqma+P0lu65a/f4AGXzNa/YZeLWOMo9my0lM1btTl2rlpvVx+LvkHBKrvZVdp6E13O5pLkrb+tVrP33GdUpOTZIxLIeERGn77QzrzvIs8nstaq28/mKyv3nlNRUWFcrn8dP6IGzTwyhvK9Bxs+2uN3n7qfm39a7WstWrVqbtGjhmnajXjPZ6rvH0+8QV9M+VNFRbmy+12q3aDRhrz2geKqlbd6Wgog7SsfN364i/6Yu5GuVxGifFRGj+6p85q6/x1jZMyC/T64j1auy9HktQsNkSjOtZS7YhAj+cqdlt9tGqfZv6VJkkK9Dca1qq6+jWKLtfMAACcSGmK2UslJVprC7wdBnDaN1Pe1MIfvtbYtz5XfINEbd+wVi/ee5OCgkN09oXDHM32yDVDFFcvQfe+/J4iq1XXqt/nafyd16lmnfrq5WC2oqIiPX79UJ194VANuf7fCggM0tyZn+vNx+9Rg2YtVa9RU4/m++nzD/TLlx/pvtdWQAifAAAgAElEQVSnqm5iU+3aslEv3zdagUFB6j9spEdzZaal6qnRw3XpjXeqx8BLVFRUqK/+86qeHj1CT06dJZefn0fzlae5M6dpxruv69ZnXlebrj2VmZaiyePG6OGRF2n8l786lou24rIb9shMJdSK1JaPrlFUWJC+nLdRl479WnNfvlRN6sY4lqug2K2Hft6u85pE6/6z6shImrk+TQ/9tE2vnJ+oIP/SNGn9zyer9uvPPTl6rn8DxYUHamNKnp6Zt1NhAS51rx/pnQcBAMBRlOYVbKUkPm5FlWet1cz3J+qGsc8pvkGiJKle4+a6Zsw4zZwyydFs2zes1d5d23XTuBcVFVtDxhi16nymLr3xTs145zVHs8368C2FR0XrH7fer+DQMPn5+6vnoMvUqc95mvrCOI/nm/n+RP3z/qdUN/FAEVw7oZGufeBpzXzf8+dg7tefqVXnHup10eXyDwhQcEioLr3xTgUEBevP3+Z4PF95+vKtl3Tx9bepbbezZYxRZEys/vXYBKXt26MNK7nyma9ZtXm/1mxN0cu39VJMRLBcLqOLejTWdQNb6Y3pfzqabcH2TNWJCNRFp8Uq0M+lAD+XBjevpgbRQZq3PdOjuYrdVl+vT9UtneMVF35gVbdRtWBd2z5OX66jVR4AULFKU8xGS1prjJlljJn+35u3gwEVrbioSKn7klWvSfMS4w2atdTeXTscSnXAlnWrFFOzlkLCwkuMJzRv5fh5oNvXr1XiaW2OaAFu1LKdUvft9ni+vbt2KKFZixJjCc1aal/SDo/PQUzetV0JzVoeMd6gaUvtTXL2Oc3NylRCs1YlxoJCQlQ9vq62rlvlUCqU1ebd6WqVGCt/v5Ivq20b19CW3RkOpTogObtQDWOCjhhvGBOs5OxCj+bKKXSr2G0V/7f25MSYICVn0cAFAKhYpSlmx0q6SNITkp477AZUKf4BAaqb2FQrF5Zs8VyxYLYaNm91jJ+qGC06dFPKniSlJCeVGF8272dFxdZwKNUB7c7spZW/z1VRYck3sktmf6d6jZof46eOLaF5Ky1fMLvE2PIFs9WgWSuPz5lteFprrVgwu0QRXFRYqFWL5irB4ec0pkYtLZv7U4mx9P17lbxzm1p1ObPCcgxNf7bEDWXTtlEN/b5mjzKy80uMf7doq9o3relQqgMaxQRr2e5suQ/7d2Ct1bKkbDWKCfZorrBAl8ID/bRuX26J8aVJ2Uqs5tlcAACcrNIUs9skLbTWzrbWzpb0u6Ty2a4UqGQuGXW7Xh97uxbMmq69u3ZozoxP9c4zD+mi629zNFdsXLyateuocTcM05+/zdGeHVs1/e1X9e3UtzT8jocczdal7wUKDg3T0zdfpQ0rl2rnpvV6+8n7tXHVMl1x+4Mezzfkhn/rrSfu09yvP9feXTs0f9aXmvToXRpyw7/LkG2gUvfu1ltP3qcdG9dp46plmnD3Daqb2FSNWrbzeL7yNOLOh/XjZ+9r2psvas/2LVq5cK4ev2GYGrVsp7g6DRzNBs/Vqxmhy85uokH3TdevK3Zq/Y5UPTh5vr5fvFXXDXT2g5N28WEK9ndp/IIkbUrN0+bUPL3wW5JcRmofH+bRXC5jdHnr6np2/k4t2J6pPVkF+m5Dmt5dvleXtWTjMgBAxTInatszxiyW1O2/G0AZYwIlzbPWdqyAfEdIbNHGjpsy04k/GqeI5fN/0Yx331DS1o2q16iZBo38l047o6vTseR2uzXp0bu09NcfVZCfp+jqNTXijod1eg/nNxbPy8nS+Dtv0KbVK+QuLlategka9dgLHm/+9F+rFs3X9Ldf0c7N61U7obEuuPpGte7co0xzZaalatqkCVoy53sFBASq27mDNfDKUQoMcn4VacX82Xrn/8Yqde9uBQQFq23Xnhr16Hi5XJ5tyHMyWI0tP8XFbr0ybbne/ma10rPzdU6H+npgRCfVj/PupkgzJs464TG5hW59unq/5m/PkLVSt3oRuqRlrEIDyrYJ2sIdmfpybYr2ZBeqYXSQLm1ZXc2qh5RpLgAA/m7w1LVLrLUdTnRcaYrZZdbadn8bW26tbXuSGcuEYhYAyo7iteorTXELAEBlVtpitjQf/+81xgz67x1jzGBJ+04mHAAAAAAAJ6M015kdJWmKMeblg/d3SBrhvUgAAAAAABzfCYtZa+1GSV2MMeEH7zt7HRAAgEdoLQYAAFXRMduMjTHDjTGHvm+tzTq8kDXGNDLGVNz1IwAAAAAAOOh4K7OxkpYaY5ZIWiJpr6RgSY0l9dSB82bv9XpCAAAAAAD+5pjFrLV2wsHzZHtL6i6pjaRcSWskjbDWbquYiAAAT9BWDAAATgXHPWfWWlss6fuDNwA4Kmutls//RXNnfq7C/Hyd3qOPzjzvIvkHBJZpvn1JO/X9J+9q56a/FN+gkfoNvUo1atcr39AAymzJrkxNXblfmfluVQ/107Xt49QwxvlrN1dmbrdb76/YpwXbM+WW1C4uVP9sH6dA/4q7rjQAVDX8BgVw0j559f/03nOP6LTTO6tTn3P164xP9dzt16q4qMjjubatX6MHRgyUu7hIPQcPlcvl0oMjLtDmNX96IXnVMDT92RI3wJtmrEvR03N3qV2tMI08vYbqRAbpnh+26s892U5Hq9Tu/G6r5m3P1MUtYjWsVXWt25+n0TM3y+12Ox0NAHxWaS7NAwDHtHfXdv3w6Xt6btocRUTHSJK69L1AY0deqD/mfK+Ovc/1aL6PXnpaF193q/oNvVqS1LHXAMXVS9CHLz2lMa9OKe/4ADzgdrv10ar9uq1rvLrVi5QkdakboWoh/npjSbJePq+hwwkrp4U7MrU7q1BvDm6k0AA/SdKZ9SP1r6836Yu1qbq4RazDCQHAN7EyC+CkrF68QG26nX2okJUkP39/dR9woVYunOvxfCsXzdWZ5w8pMXbm+Rdr5e9zZa096bwAym5vTpFyC93qUjeixHivhEjtzS5wKFXl9/OWdHWtF3GokJWkAD+jPg0j9duOTAeTAYBvO+HKrDEmSNIQSQmHH2+tfdR7sQD4irDIKKXt3XPEeOrePQqLjPJ8vogope7do9Dw/71ZTtuXrLDIKBljTiorcCoYeH3/EvdnTJxVbnOHBbpkJWUVFCsy6H9vIVJyixTg4vPxY4kK9FNSVuER4/tyihQe6HeUnwAAlEZpXnm+lDRYUpGk7MNuAKC23c7W7m2btWDW9EMrp1vWrtScrz5Rjwsu9Xi+sy8cpg/GP6783FxJUkFert5//lH1GjysXHMD8Fx4oL/qRAZq0h/JKiw+8O89u6BYk5cmq2XNEIfTVV7DWtfQmn25WrQz69DY2n25mrM1Q8NaV3cwGQD4ttKcM1vXWjvA60kA+KSAwCDd+cLbmnD3KH0x+SWFhIUraesmjRwzTvH1PT9/7uLrbtWkx+7Rzed3VkKzVtq6bpVad+mhIaNu90J6AJ4a27Ou7v1xm66atl51o4K0JTVPdaICdUfX2k5Hq7RiQvx1bfs4PTt/p2KC/RXgZ7Q7q1AXNa+mprF8CAAAZWVOdA6aMWaipJestZViK9HEFm3suCkznY4B4G/cxcXasHKpCvPz1aRNewUGn9wbtL27dmjXlg2Kb5ComnXql1PKqoEdi+GJ8mwzPtwfu7K0MTVfp9cKVWMKslIpKHLrh81pKii26pcYpdBA9uEEgKMZPHXtEmtthxMdd8zfosaYPyXZg8eMNMZskpQvyUiy1to25RUWgO9z+fmpadsT/s4ptRq166pG7brlNh+A8tW+drja1w53OoZPCfR36bwm1ZyOAQBVxvE+EhxYYSkAAAAAAPDAMYtZa+1WSTLGvGetHXH494wx70kacdQfBACUK9qKAQAAjlSa3YxbHn7HGOMn6QzvxAEAAAAA4MSOWcwaY8YYYzIltTHGZBy8ZUpK1oHL9QAAAAAA4IjjtRk/KelJY8yT1toxFZgJAE5ptBUDAACc2PF2M25/8MtPDvv6EGvtH15LhVNCSnKSNq/5U9VqxiuheSsZY5yOdMjOTeuVtG2T6iY2Va0yXCv1cPv3JOnnaR8oJDxCfS+5UoHBweWUsnIpLirS2j8WqqAgX6e176zg0DCnIwGlsj89V/NXJalaRLC6toyXy1V5fhedKvZkFejHTemKDPJTv0bRCvQvzVlQR+e2Vuv25SqzoFjNq4coMqjqXv5me3q+dmUWqH5UkOIjAk9qrpzCYq1OzlWQv1GLGqHyO4l/B9Za/bU/T2l5RWoWG6LokMrzHBS7rVYl56jAbdWyRqhCAsr+/xoA5x3vt8tzB/8bLKmDpOU6cFmeNpIWS+rq3Wioqtxut6aMf0xzvvpETdqcoZ2bN6hajTj9+7lJioyJdTRbXm6OXrn/Zm1cuUwNT2ujjauWqVXnM3XD2GcVEBjk8XyTHrtbc7/+XA2atlBOVqY+e/15jRwzTj3OH+KF9M7ZsHKpJtw9SlHVqis4NFyvPnCLrrr7MZ153kVORwOO6/mP/9C4935Xx+Zx2rUvW9Zaff74QDWpG+N0tFPGs/N3aeGOTDWqFqy0vGJN+XOfbusSr851IzyeKymzQE/+ulNua1U9NEAvLEjSJS1idXELZ19byltekVvPzd+l9Sl5ahwTpL/256ldrTDd3DleAX6eF6E/bErT20uTlRgTrOxCtzLzi3XvmXXUqJrnH77uzS7U0wuS5PbzU0KtSL28eKvObxKtoS2qOf6h9bp9uXpm3k7FBPsr2N+lFxbs0vVn1NJZCZGO5gJQdsdrM+4lScaYzyW1t9b+efB+K0kPV0g6VEmzp3+sdcsW6YWv5iksIkput1tTJzyhSY/epTvGv+Voto9eekpBwSF68evf5B8QoIL8PE24+0Z9MfklXXrjnR7N9fuPM7Xgu6/0xAffqE5ik0Njrz14mzr26q/g0KpxfcaC/Dw9f/u1GjlmnDr2GiBJ2rFxnR6/fqgaNm916LHj2GgrdsaPS7bplWnLtWzyFapXM0LWWr32xQpdOvZrLX3zCsffeJ8KvtuQpqVJ2XrpvIaqFX5gZfGHTWka/1uS3r0wzKMVWmutnpm3U30bRWlg0xgZY7Q/p1D3/bhNCTFBah9fNX7nStI7y5IV7O/Sm4Mayd9llF/k1tPzdurT1ft0eesaHs21OTVP7y3fq6f7NlDdyAMf2s7blqFxv+7QGwMbeVwcT1i0RyMuaKP7hneSMUbJqTk6a/RHarAjS13ref4BRXnJL3LriV936KZOtdSpzoEcW9Py9cBP25RYLejQYwfgW0rzKtHsv4WsJFlrV0o6zXuRUNXNnv6xLhl1h8IioiRJLpdLl4y6Q2uWLFRG6n7HcllrNeerT/SP2+6Xf0CAJCkwKFiX3zpGs6d/7PF8M959Q/2HXl2imOvU5zzVTmisr955vdxyO235/F9Up2GTQ4WsJNVt1Ew9Bw/Vr19/5lww4ATembVGd1zWXvVqHnhja4zRjRe2UV5Bsf74K9nhdKeGmRtSdfFp1Q4VspJ0TmK0YoL9NWtjmkdzbUnLV06hW+cfLGQlKTY0QBefFqufNqWXa24nua3Vz5szNPL0mvI/2Aoc5O/SVW1r6scyPM6fN6erf+PoEsVc9/qRqhkaoOV7sj2aKymzQHtyinTPPzoeeg5qxoTqwau7avaOLI+zlaclSVlqEBV0qJCVpAbRQerTMEq/bMlwMBmAk1GaYnaFMeZNY8zZB2+TJK3wdjBUXblZGUe0EwcGBysoJER5OZ69cJYnd3Gx8vNyFRFdrcR4ZEx15WZlejxffm6OoqvXPGI8unpNZaamlDlnZZOblanIake28EVWi1VutrNvXoDjycguUI2Y0BJjxhjVjA5RenaBQ6lOLYVuKTrY74jx6GB/pecXeTRXTqFbkUF+cv1tRT0q2E85he6TylmZuK1UUOxWRGDJv7eyPs6cIreijnJecVnmyyl0KyY8SP5+Jd9e1owJcfw5yCl0Kyq4fB4ngMqjNMXsSEmrJN168Lb64BhQJq27nKU5X31SYuzP3+YoODRM1ePrOpRK8vP312lndNWvM0quJs6e/rFad+3p8XxtuvbUj59NUXHR/96Qpe7do1WL5qn3xVecdN7KomXH7lqxYLbSU/YdGisqLNS8mV+oTdezHExWeQ1Nf7bEDc7o26G+3p21WtbaQ2Nrt6VozbZUdT6tloPJytfA6/uXuFUmzWKD9N3GdLkPew72ZBVoQ0qu+iZGezRX42rBSsos1Lb0/ENj1lr9tCldp8dXnQ3p/F0HNmiavaXkKuyPm9LVrgyPs12tMP2yJV3F7v89B/tzCvXnnhy1rhl6nJ88Uv2oIO1Lz9XidXsOjVlrNXnGSrWu7uzmh23iwrQ0KUvpef97TS4stpq9JUOn16o6/38Apxpz+Iu4L0hs0caOmzLT6Rg4Cekp+/TINRcrsUUbtT/rHO3YtF4/fTZF/3p8gtqUoWgsT1vWrdJTNw1X9wGD1bh1e61eskBLfvlOD0z8WLUTGnk0V0Fenu64qKciYqqp/7CRys5I11f/eVUNT2utu19610uPwBmfT3pBc776RAMu/6eCQ8P087SpiqwWq3//30S5/I5cdTnVUcBWDjl5hep35zRFhAbqir7NtWtfll78bJkevaarrjmvpdPxvGbGxFlORzgkp6BI/5q5RXFhAerfOFrpeUX6bHWK2tQK1d3d63g830+b0/XusmQNbFpN1UP9NWdrhjIKivV47/oKPokdkiubjSl5euSX7eqZEKmmsSH6MzlHv+/M0hN96qu2h7saF7utxs3Zodwit85JjFJ2gVtf/ZWiAY1jNKQMG2fN25ahySv2a/SQdmpcJ1of/rBWazfu1WM96ygs0NnXgw//3Keft6RrYNMYhfi79N3GNMWE+OueM+scsaIPwFmDp65dYq3tcKLjjlnMGmM+ttZeZoz5U9IRB1lr25x8TM9RzFYN2Znp+uWLD7Vh5TJVq1lLvS++QnUaNnY6liRp/+5d+vGz97Vr6ybVa9RMfYZccdR24dLIy83RlOcf06pF8+QfEKizLxym8664tpwTVw6rfp+ned9MU0F+nk7v0Udd+l4gP//KczmGyoRitvLIzS/S+9+v1Y9LtqlaZLBGnttCHZtXnVXZo6lMxax0oKB9c2my1u7LU4DL6NzG0RrQpOy7SW9IydMPG9OUUVCs1jVD1bthlIKqUCH7X3uzC/XthjTtyixQg+gg9W8UrZgyXgKnyG01d2uGFu/KUpC/S70aRqmVh6uyh9ucmqeftmYqLd+t5jGB6t0wqtJcAmf57mzN3pKhQrdVx9rh6l4/4qQuQwTAO8qjmI231iYZYxoc7fvW2q0nmbFMKGYB+BqKV1Qmla2YBQDg70pbzB7v0jxJB788R9Ica+368goHAAAAAMDJKE0/Sn1JbxhjEiQtkTRH0q/W2mVezAUAAAAAwDGd8AQGa+1Ya21vSS0l/SrpLh0oagEAAAAAcMQJV2aNMQ9I6i4pXNJSSXfqQFELAAAAAIAjStNmfLGkIklfS5otaYG1Nv/4PwIApzY2fQIAAPCu0rQZt9eBTaB+l9RX0p/GmLneDgYAAAAAwLGUps24laQeknpK6iBpu2gzRhW2ZPZ3+uqd15W0ZaPqNW6mwdeMVusuZzkdS9Za/fDpe/rhk/eUvn+vmrXrqCGjblf9JqeVab75s6brw5eeVMb+fQoMDlHbrj114+MT5HJVjmsBAjg17Mkq0Ad/7tOy3dkKDXCpT8NoXXhaNflXsWt/JmUW6MGftyszv1jWWoUH+eue7vFqVr3s13MtL7/vzNQrv+9WbqFbLiNFB/vr6XPqKyokwOloAHBcpWkzfkoHdjB+UdIia22hdyMBzvnt+xmaMv4xXXX3o2rSur3WLPlNrz30b90w9lm17d7L0WyfvPasVsz/RSPvfVxx9Rpo4Q8z9cSoyzX27WmKr9/Qo7mW/vqTJj5ypy6/ZYw69TlPu7dt0uRx9+nJG/+h+9/40DsPoIqjrRjwXHpekcb8uE39EqM1vF8NpecX673lydqZma9bu9R2Ol65KS4u1h3fbVXbWqG6vFV1BbiMvliXorE/79AbAxs6WjRuTs3Vs/N26fymMRrQOFrZhW79Z1myRn+zWe9d3NSxXABQGqVpMx5orX3GWjufQhZV3WdvjNeoR8arw9n9FRVbQ136XaCR9z6uaW++6GiunMwMfffRO7pj/Ftq3r6zYmrU0oDLr9E5l16pme9P9Hi+D196QoNH3qT+w0YqpkacTjujq8a89oH+WrFYKXv3eOERAMCRvtuYptNrhWlY6+qqERagxtWCdV+Pulq8K1tJmQVOxys3k/5IVniAS3d1q636UUGKjwjUqDPi1CAmSM8tSHI02wsLktSuVpiualdTceGBSowJ1oNn1ZPbStPXpjiaDQBOhH5C4KCiwkLt2rJBLTp0LTHeqnMPbf1rtUOpDti9fYtq1K6rmBpxJcZbdT5T29ev9Xi+9JT9atW5R4mx2Lh4RcfW1Ppli08qKwCU1ua0fLWNCysxFuTvUrPqIdqWXnX2mly9N1en1w6Ty/yvddoYo461w7XT4aI9vaBYHeqElxgL8DNqUSNUi3dlOZQKAEqnNG3GwCnBz99fsTXjtXXdKiU0b3VofNPqFYqr28DBZFL1+Drau2uHcjIzFBoReWh88+oVqlmnvsfzhYSFa9Pq5WrSpv2hsayMNKXt36uE5i3LJXNVR1sxfNXA6/sf+nrGxFkOJpHiwwO1MTVPZyX87/dasdtqc2qe4sKrO5isfNWPCtLavblHjK/dl6uYYGffioX4ubRuX676NYo+NOa2VhtS8tTrsOcFACojVmaBg4wxOv/KG/TGw3do+4YDq52bVq/QW0+M0fkjbnA0W2RMrDr1OVevPXSb9u9Jktvt1tJff9SXb7+iAf+4xuP5Bo28SR+/8oyWzftZ1lrtS9qpCXePUnyDRMXVSyj/BwAAR9GvUZR+3pyuXzanq9htlZFfpFcX7Va9yCAlRAc7Ha/c3NSxlvZkFeq9FXuVU1isgmK3vlybohW7s3Vbl3hHs13fIU5ztmbouw2pKiy2yswv1uuLdiuvyK0rWsc6mg0ATuSYHwcaY76SZI/1fWvtIK8kAhzUb+jVchcX66mbhisnK1MRUTG68Npb1GPgEKejaeS9j+vjV57RPZf1VVFhvmrVa6gbHx2vxBZtPZ6r14XDlLInSa/cf7PycrLlcrmU0Ky1Hn7rMy8kB4CjiwsP1AM96+qtP5L1yqLdchnpzPqRuvvMqrP5kySFBPrprm7xmvD7Hk1bvV+SFBbop+va11TdqCBHs50eH65hrWL1zvK9en3xgT0TwoP89FivevLz83M0GwCciLH26PWqMabn8X7QWjvbK4lOILFFGztuykwn/micQtxut/JyshUcGlbpLlVTXFSkgvw8BYeGyZiTu3SF2+1W2v69Co+IUmBw1VkF8QbailEVOd1mfLi8Irf8jFGAX9W6JM/fpeQWyFopNjTQ6ShH2JNVoBB/P0UGU8QCcNbgqWuXWGs7nOi4Y67MOlWsApWBy+VSaHiE0zGOys/fXyH+4Sc+sBRcLpeq/W1TKQBwQrB/5frg0FuqhVS+Iva/4sIrbzYAOJoT7jpgjGki6UlJLSQdWrqx1iZ6MRcAAAAAAMdUmi303pY0VtJ4Sb0kjRQbRwGo4mgrBgAAqNxKU5SGWGt/1IHza7daax+WdL53YwEAAAAAcGylWZnNN8a4JK03xoyWtFNS+ZywBwAAAABAGZSmmL1VUqikW/T/7N13eJRV2sfx7zOT3hMCJARC70iNCEoTFESxF2xrF+vay7s27Gtbd+0rrg3FiopSpErvHeklAZIQQnovU877RzAaQWVCwqT8PtfFJTnO3NyTZ9r9nPs5B54BhgPX1mZSIiInmtqKRUREROqXvyxmjTGrAQ7Pzt5ljCk4lsCWZQUAiwD/w//OZGPM+N/dxh+YCPQDsoCxxpi9njwAOXFyMw+xdOYUCnNz6JYwkO79B1V72xqX08m6xXPZs3kDUc1iOHX0BYSERdRwxg1PatIuvnzzJfKyM+g5YCjn3/h3fHyO5ZzUkdwuFxuWzmfnxjVERDfj1LMuICwyqtq5pafsY/msH3CUl9Fn8Ag69OhT7ViNhdttmLNmP4s2ptA0IpArRnSmeVSwt9MCIL+wjEffX8b6XYfo2DKSf447jZjjyG3j7gymLNmD3WZx8dCOdG1d/edaTXK73fz76/VMWbybqLAAxl83gL6dqr/Cd0pGAZ/P3UFuYRlnntyaob3iqr2FVnFpOY+/v4IVW9Jo2yKcF245jZZN6+Yq63XJ+rRC3l93iBKnm97Ng7nl5Gb4VXO/VJfbsCq1kF3ZpUQH+TCkdRghftXftmbtgQKmbM/B5Tac1SGCIW3Cqx2r1Olmyb58DhSUEx/hz6mtQvGzN7wlTVxuw9q0QnZklhIZaGdI63DC/Kt/DNIKylmyP59yl+HkuBA6NQmswWwbJrcxbEgrYktGCeH+doa0CSMioHrfPaBiC6jF+wsoc7rpGxtMl+jA495qUBqvv3zXsywrwbKsn4FNwM+WZW20LKvfMcQuA4YbY3oBvYGzLMsa8Lvb3AjkGGM6ULHA1IuepS8nyuaVS3josjNI3bMTH18/Pn31af79wM04HQ6PY5UWF/HMuMuY+uHb+AcGsnPjGh68eDh7t2+uhcwbjjlfT+TRq87B18+fXgOHsWzmFO4eM5DiwnyPY5WXlfLPO65m8n9fxS8ggMStm3jokuHs2rSuWrktmvo1j19zLnlZGbhdLl5/+HYmvvwkf7SPtUC5w8UFj03l/95dgr+fnZ8Ts+h5wyQWbEjxdmps25dF2ys+ZN3OQ5w9oC0Hs4vodNXHzF+3v1rxnv54BWP+8T0lZU5yC8sYfu9k3vhmQw1n7bnycicdrviQVz5ZSSvLQUl6DoNu/5JnP8XCWwEAACAASURBVFlZrXhTlyXS87pPWLhgG0mbkrj+mRlc9fQM3G7PXwf7DubT8qL/MW3+Ntr7OtmzI5XOV37EtOVJ1cqtsXhjRRrPL06lQ1QAZ7QLZ9OhYm6YsoeScpfHsUocbh5fmMrc9DJ6DehEVmAId83aR2JOabVye3lpKi8vTSM+wp9O0YFMWHuIR+dV7zWVXljO32cksSK1kABfGz8l5nHvzL3kljqrFa+uKnO6eXJBMl9szsLPbrErq5Q7pyeyI7OkWvHmJuby0Jx95Ja6cBt4aUkq769L12fVn3C4DM8uSmHixgx87RZJuWX8fUYSWw4VVyveon353D9rH1nFFd8f/7MijXdW6xhI9Vl/9eSxLGsTcIcxZvHhnwcBbxtjeh7zP2JZQcAS4DZjzMrfjM8CnjTGLLcsywc4CDQ1f5JUu249zXOTZhzrPy01wOV0cve5p3LbU/+me//TAHA6HDx/2xUMHnMJp19wuUfxvp3wH5L37ODv/3yrcmZ38bRvmPXFBzz76fQaz78hcDqd3DK8J7c/8xr9hp4JVMysPnfr5YSER3LvKxM8ijf9k3fZuno59//7fWyHZyxWzfuRr99+mZcmz/PoDGlBbg73nj+Ypz/+nhZt2gNQXJDPo1efw7gnXqJrv4Ee5eYtJ7rN+L8/bOKbBbv48eUL8Tk8mzJr1T5ue3UeuyZdh92LMyx9bprEaT1iefOe4ZVjL362mv/+8DNJX9zgUaxNezI4++HvWf+/K2kaEQTA/vR8+o37nDXvXkHrmLAazd0Tt/5rHvOW7ualM+PxPfz73nKomKcWJnPo+1sJCTr2PTdLy53EX/w/Hh4YQ+foipmecpebxxak8tRtQ7l0WEePcjvlls8JdZRx9ykxla/H6btymLIzl7Tvb/Uo1p+ZNmFWjcXytsJyFzdM2c3jQ1tyUvOKLgKHy/DwnH34+cALZ7TxKN4XmzNxRIbzxZPnVB6Dj2du5eWPlvHC8FYexdqZVcJj8/bz5jntaBbsC0B+mYvbpyVyc79mDPVwhvb5xSl0jArg0u7RlWPvr0un1Onmjv6xHsWqy77dlsXWjBL+MSgOu63iGCzdn88XmzN5fXRbjz6r8kqd3DYtkZdHtiEurOK1XVTu4r5Ze7nrlFi6NwuqlcdQ383YlcOK5ALGD2tVeQzWHijk3TXp/Pfcdtg8OAZF5S7GTd3DcyPiaRNRsdtnicPNg7P3ckOfZvRtoSV55Ffnf759rTEm4a9udyzflly/FLIAxpglwDGd+rMsy25Z1gbgEDDnt4XsYXFA8uG4TiAPaHIsseXESdr2MyFhEZWFLICPry8jx17Hmp9mehxv9fyZjL7yxiotyqeNvoDMtFSyDh6okZwbmvWL5uAfEEjfIWdUjtnsdsZcexu7N6/3ON7q+bM468obKgtZgJOHn0VZaQkH9u7xKNbGZQvonjCwspAFCAoNY9j5Y1n9U8P5olzTvl+SyB0X9a4sZAFG9W9NUIAvG3ZneDEz2J2ay/1jqzbg/P2i3hzILCS/sMyjWFOW7OGqM7pUFrIA8c3DuGBQe6YuS6yRfKtr9sokzu8SWVnIAnRvFkSzYD8+nLnVo1iLN6XSIsyvspAF8LPbGNkmlK9/2uFxbtv2ZnFBl8gqX9ZHtosgp7CcfQc978ZoDCZvySQq0KeykAXwtVtc2DWK1HzPu4hWHyzh3rH9qhyDq8/swsFCR+Ws0rH6YXs2p7QKrSxkAcL87YzsEMHMXbkexXK5DWsPFDGmU9VW/fM6R7E8udCjWHXdipRCzu0UWVlEAZzaKpRih5sDBZ4dg3VpRfSKCa4sZAGC/eyMaBfOipRjuoKuUVqRXMA5vzsGfWOD8bFZ7M3x7PNgU3oxnZoEVhayAIG+Nka2j2BFSsN67sqJcywN7wsty3oX+BwwwFhggWVZfQGMMX/Yl2iMcQG9LcuKAL6zLKuHMcbjXlLLssYB4wCiY+I8vbscJ5vdhsvpxBhT5UPd5XRg8/H8uhWbrSLebxljcLtd2Brg9T41we7je8TvDCqOgYXn15nYbDacRzsGLid2D68ts9ltuFxHy81Zp4+ntxd8stssHE53lTFjDA6nq0qB6w0268jcnK6Kn202z55vdpsNh+vIFk+H0+3V2WcAy7JwHaUF2Ok2+Pl4lpvdZsN1lKYipzH42j1/jVoWuKseAtzGYAB7NeL9kTHjRlX5uT7P1PrYraMeA5fbVONdEmzWr8/7X7iNwW0Mdg+v77PbLMr/4Llm9/A19Utuv3+sLmOow2+51WI/yuN0m4pj6uljtVkc9fXuchuPZhcbG5tlcbQrJVzG4OmyKX90DJzGUI2XgQhwbDOzvYBOwHjgSaAr0Af4F3BM3waNMbnAfOCs3/2vVKAVwOE243AqFoL6/f0nGGMSjDEJocexQI1UT5suJ+F0Olj9m1nY0pJifpz0PwaOPNfjeANHnccPH71d5XrbuZM/oWW7TkQ2jamRnBua3oOG43a7WDz9m8qx8rJSvvvf63RL8LyNd+Co85g+8b84yn89q7roh68Ii2pK81ZtPMvttOHs2LCaxK0bK8fysjL46bvPGFCN50djcdnpnfj3V+soLv31dfDV/F342G30bB/9J/esfV1bR/HUxyuqXOv5/Kerad08zKPWW4BLh3Xks7k7qswmbtuXzbTliVwwqP2f3LP2XTSsI99sy6bY8WuxvTq1kNxSJ9ef1c2jWIN7tiCrxMWGtKLKsaJyFz8m5nPFyK4e59anU3O+2JJZ5Yvf9zuyaRYeoEWg/sDFXZuQX+qqMstW4nDz9dYs2kcF/Mk9j25giyCen7gSh/PX58db322kTUQAEYGeLX5zcbco1qQWsj/v1/fcjCIHs/fkcm7nSI9i2W0WA1qGMnlLVuV1hsYYvt6SxaB477Xt14ZB8WF8uzUbx29OKsxLyqNJkC8xIZ69F/VrEcKWjGL2ZP96zXNOiZM5e/IYFK/X1B8Z3DqU77ZnU/abE5yL9xXga7NoHe7vUazeMcEk5ZaxLePX623zy1zM2p3LoNYN67krJ85fXjNb7cCW1RRwGGNyLcsKBGYDLxpjpv3mNncAJxljbrUs63LgImPMZX8WV9fMesfuzev517030q5bL6Jj41i7cDY9Bw7lpsde9HhFY0d5Ga89dCupSbvpfdrppCTuIj15L/945zNi49vW0iOo/1bMmcq74++nXfdetGjbgdXzfiQoNJyXvpqDj59nH+oup5O3Hv07e7Zuos+g4Rzcn0TKnh08/OYntOrQxePc1iyYxbtP3k/PgUMJCAphzeFW8gtuusvjWCeKt2dm3W7DjS/NYf76ZM49tR17D+azbuchfnj+PPp1rv5qujUhJaOAk8d9TmiQH6P6t2bpzwfYl17AvH9fTO8OTT2O99Z3Gxn/4XLOP609Dpeb6cuTeO2uYVx9pufPtZrkdrvpe+Mkkg7kcVp8KBnFTjanF/Gfu4Zxy3nHvCxEpUUbU7n4sR/o1jSQcD87K1MLGTuiM6/dfbrHK3Vm5hZz0rUTwe0moUUIu7JKSS0oZ9YrF3HqSS08zu1Y1eeZWYBPNh7ihx05dI0OpHmIL8uSK750TzivnccrGjtcbl5ecZBDZYazB7Rh055MdidnM35wHLGhnr3nAry75iBzE/M4OS4EX7uN5fvz6dcihIcHed5xllvi5In5yQT4WHRqEsjmQ8X4+dgYP7Qlwcex2nJd43Qb/rXsAHuyS0mICyE1v5zk/DKeHNaKeA8LKYAVKQW8sTKNPjHBBPjaWJFcyAVdoriku65w+yMut+H1lWlsOVTMyXEhpBc5SMwp44khLWlXjZNE6w4U8uryA/RsHkyIn50VKQWM7BDB1T09/2yRhu1Yr5k9lgWgmgPPAy2MMaMty+oGDDTGvP8X9+sJfAzYqZgB/soY87RlWU8Da4wxPxzevucTKmZ6s4HLjTF/ehGVilnvKS0uYvVPMynMy6FrwkDadO5e7VjGGHZuXFOxNU/zWPoNPRNfP88/mBqb/Jxsvnn3X2SnHyRh+CiGnvun537+lDGGPZs3VGzN07QZCUNH4hdQ/S0K8nOyWTN/Jo7yMnqfdrrHM7y1zdvF6x9Zt/PQ4a15grhgUHuCA33/+k4ngNPp5tWv17Fs8wF6tI3msWv6E+BX/a0Y9qfnM21ZEja7xQWD2h/XNj81bfKCnXw2bwdNw4MYf90ptIiu/iIkeYVlfLdkD/mFZYzoF0/3ttX/kux2u3n9240s3JBC51aRPHFtf4ICPC+iPFHfi1mo2HrlrVVp5JW5GNYmjIu7Vb/TwRjDtswSdmVVbM3T/3AhWl37ckqZvC0LpxvO6xxJ16bVX3Toly1r0gocxIf70ysmqEG2yxpj2JlVyvbMEiIDfTglLgR/Dy8D+K38Micrkgspd7vpFxtSrRMTjdHu7FK2HComPMDOgJahBBzHMSgoq+igKHW66RsbUuU6ZpFf1GQx+yPwIfCoMabX4Xbg9caYk2omVc+omBWR6qirxaxIXdIQilkREan/anI142hjzFeAGypXHfZ8wzYRERERERGRGnIsPWNFlmU1oWIlYyzLGkDFFjoiInWaZmNFREREGq5jKWbvA34A2luWtRRoClxSq1mJiIiIiIiI/Im/LGaNMessyxoKdAYsYIcxxvPdx0VERERERERqyF8Ws5ZlXQrMNMZssSzrMaCvZVnPGmPW1X56IiLHTm3FIiIiIo3HsSwA9bgxpsCyrEHACOB94J3aTUsaA6fDQXryXooL8r2diogcRWm5k92puRSV1EwzTkpGAQcyC2skVl2WmVdCUloebvfx7+NedvgYFBSX10BmkJpRSGpGzRyDEoebtIJyHC73ccdyG0N6YTl5pc4ayKxmOVyGtIJyih1a+1KOVOaseB2UOY//dSAinjuWa2Z/efc+B3jPGDPdsqxnazEnaQTmf/c5X7/zCj6+fhQX5DNg5Llc8+CT+Pl7vgG3iNQsYwz/+nIdL32+hrBgP3IKyrh5TA+eu+lU7NXYY3PD7gxueWUu+9ILcLkN3VpH8d6DZ9CpVWQtZO89mXkljHtlHgvWJxMc6EuQvw+v3TWMs/q3qVa8N7/dwLOfrCI4wJecglKuGdWNl28bhK+P3eNYW/dmMe6VeexMzsGyLNq3CGfCgyPo0dbzPVidbsPHGw4xLzGPYD87ZU43l3RrwnldojyOBbAhrYj3N2VS5jIUlznp0TyYW/s2JSKg+vsa15TZe3L5cms2Af4+5BWVM7h1GNf3isbvOPaalYbBGMOXm7OYujObYF87RQ4XYzpFMrZHdIPc71ekrjqWT4pUy7LeBc4EXrQsy59jm9EVOar1i+cx5f03+L83PyW+U1cK83J475mH+eSVp7jx0X96Oz2pR9RWXDs+mLGFSXO2s/ztsbSPiyAtq4irn53Js5+sYvx1AzyKlVtYxjkPT+Gf4wZx1RmdcRvDhKmbOevB79g68RoC/LxfsNSUS5+YTr/Ozfj00ZsI9Pdh3rpkrn52Jj/9+2K6tWniUawvftrB21M2seC1S+gSH8WhnGKue2E2j/5vGS/dOtijWMWlDkY/NIVHru7Pjed0xwI+mrmV0Q9NYevH1xAa5Fd52zHjRlW579H2nf1sUwb788p4+5x2RAT6kJpfzvOLUwjztzOsbbhHuaXkl/Gf1el8Nv5szkyIp7TcxfgPlvHSol08NywOy4tFwcqUAqYmFjD3tUs5qV00WXkl3PjibD7alMW4Pk29lpfUDVN35LAmrZD/nNWWpsG+ZBQ5eGlpKoE+Ni7o6tnrXUSq71iK0suAWcAoY0wuEAU8WKtZSYM28/MPuPzv/0d8p64AhIRHcvPjL7J81g9qORapA974dgP/+ftQ2sdFABDbJJh3HxjBW99t9Lh19vN5OxjSqyXXjOqK3W7D18fOHRf2omOrSL5fklgb6XvFxt0Z7D9UwIu3DCIowBfLsjijXzy3nd+TCVM3exzvjW828Mrtg+kSXzHb2SwyiPceOIP3p2+hrNyzVtzJC3fTq0NTbjnvJHzsNux2Gzee04OTu8Tw9YJdHsVyug2z9uRyZ/9YIgIrTkTEhflxU9/mTN2Z41EsgLlJ+dxyfk9Gntway7II9PfhhVsGU+iGPTllHserSbP2FvDSHUM5qV3F7HWT8ED+9/BIFu3NU8uxMHVnDrclxNA02BeApsG+3HZyDNOq8ToQker7y2LWGFNsjPnWGLPr8M9pxpjZtZ+aNFTZhw4S165jlbGQ8EiCQsMoyM32UlYi8ovUzCK6tanaMtq+RTiFJQ5KPSykUjMK6db6yPbTbq2jSG1A18+mZBTSOT7yiDbsrq2jSMko8DheambhEbO5LaKDsSzI9/D62QOZf3IMPLx+ttTpxm0gOqjqjHqrcD+yij2/tjqnzE33tlUfp81m0bllRLXi1aSsYucRv7fo8EBCAn0pKFMx29hlFTtoFe5XZaxVmD9ZJU6MOf7r5UXk2KhdWE64Dj16s27hnCpj+3duw+kop0lMnJeykvpgbN4rVf5I7Ti5S3OmLkuqMjZnzX46xEUQ6O9ZW3D/rjFMX5FUZUbX4XQxc9Ve+ndtXiP51gV9OzVj1baDZOeXVhmftjyJU7rGeBzv5C4xTF1adeZ62eY0IkMDaBIW6GGs5vy4ci+u3yzU5HK5mbEiif4e5hbsayMiwM6WjJIq46tSC+nUxLO8ANqF+fL94t1VxnILy1i5PZ32Ud5dQ6FDpD8/LN1TZWzD7gycThdNgny9lJXUFZ2aBLI6terJoNUHCujYJNCr7fEijY39ySef9HYOHnntrXeeHHHxVd5OQ45DXNuOvP/8P3CWlxMUEsaWNct475mHuHjcvbTv3svb6Ukd1qNsmbdTaBQ6xEVw00tz8fexERzoy/QVSdz9+gJevWNoZdvrsccK57O52/lx1V7iokNIPJDH319fQHR4EA9dkdBgvvSFBvmRmVfCsxNXEdc0hIJiBy9MWs2C9Sm8c99wj08CdG4Vyc2vzMNmVcSevXofd/xnPs+PO62y7fVYtY0N4/uliXyzaDctm4aQfKiA+95ajN1mMf66AX96DHaurVrMWZZFZIAPb606SJifHcuChXvz+WpLFnf2jyUq0LPH2SrMj49WprAzNY+YqGA27snk+udn0auJP6e1CvUoVk2LC/XluWnbcbrcRIQG8NO6ZG56YRaXdo6kg5cLbfG+ZiG+vLHyIAE+NvzsNpYnF/DR+gxuSWhOTIjfXwcQkT/1xebMtCeffHLCX93Oqm+tEO269TTPTZrh7TTkOB3Yu4epH73N7s0biGoWw6jLr6fvkDO8nZbUcZqNPXHW7zrEy1+sZdOeTDrERXDvpX0Y2rtltWIVlTh49at1TFmyB7vN4rLTO/H3i3rh34AWf4KK1U0//HErH8/cSm5hGWf0i+fBK/oRExVcrXibkzJ56fO1rN95iNYxYdx9SW/OTGhdrVglZU5e/2Y9kxfsxmC4eEgH7r6kD0EBfz7DeLQFoAA2Hixi6o5sDhY6aBsZwMVdo2gTWb0CL7fUyZQdOWzKKCHI186QlsGc0S68TqwIm5JfxpSduezOKSM6yIez2oaR0CLE22lJHbEjs4TvtmeTnFdGqzA/LujahC7RnncoiMiRzv98+1pjTMJf3U7FrIjUGypmRU6sPypmRUREatOxFrO6ZlZERERERETqHRWzIiIiIiIiUu80rAuWRKRBUVuxiIiIiPwRzcyKiIiIiIhIvaOZWRGpMzQTKyIiIiLHSjOzIiIiIiIiUu9oZlZERI6wJzWXf3+1jo27DtGxVST3XNaXnu2bVitWYXE5t/xrHks3H8CyLEb0i+ftu4fhV819ZldtO8gbk9ezNy2PU7rHcvelfWnVLLRaseqyZT8f4P63F5GaWUhYsD//d2UCV4/sWq1YTpebj2du5cu52zHGcNmILlx/dnd87NU7p52YXcq0nTmkF5XTNjKAcztF0jzEr1qxsoodvLnqIPtyy/CxWQxuHcrfejWrViyAdWmFzN9fSEG5m25RfpzdMZIQP3u149VV+3LLmLEnj4PFTlqF+HJOh3BiQ6t3DMT7yl1uZu/JZVVKIX52iyGtwxncOhSrDuy3XJOMMSxPKWB+Uj6lTjf9YoM5q2MkAT6aX5Pq0T6zIuI1aiuum35OzGTEPZMZ0TqUk5oGsjunlKm78vjqmTGc3qeVR7GcTjcdr/qIuOhgHrgiAafLzfOfrKK03MnWidd6nNuUJXsY9+JsLugUSetwP9anF7M0tYilb19OuxbhHserq2at2scl46dx63k9GX1KG9bvPMQzn6zkvkv78sR1AzyKZYxh7Php7Nh9iDHtK35HMxLzaNM6mm+fO8+jL8vTJsxiXVoh/1mexgVdo+gQGcD6g0X8lJTHcyPiaRnm71FuuaVO7pieSLemQZzVIYK8MheTNmXQItSPZ4bHexQLYNrOHGbtK+Sx6wbQqlkon87aypL1yTx/essGVdBuOVTMyysOcv/YfpzSPZb565J5Z8oGxg+Oo21kgLfTEw+53Ibx85PxtVuM7hhBqdMwZVsWnaMDuSUhxtvp1aiJGw6x+kAhF3dtQoi/ndl7csktcfLsiHj8qnlyTRqmY91nVjOzIiJSxePvLeHCjuGc2zkKgJ4xwcSG+PHAmwtZ+/7VHsX611dr8bFbzH/tEnx9KoqJcwa0pe3lH/DJ7G38zYOZRrfbcN8bC7i3f3NOah4MQK+YYPztFs9+vIIP/jHKo9zqsnvfXMj/XZnAo387BYDhfVvRo100Vz7zI49d0x+b7di/9C3fksbKzQf4z5nx+B7+spjQIoT75+5n8aYDDOkVd8yxjDF8tD6Du06JJSEuBKh4foT62fni50weOO3YYwH8b2067SMDeGRwXGVR3TcmmJun7mFfTimtPSjMih0uvtyazbr3r6JtbEXRPvLk1lz51Axm7c7l4m5NPMqtLpu0NZu37x/BZad3AiqeH7HRwXw2ZT3/d2qsl7MTT61MLaDc5ebp4a2xHX4dJLQI5tapiZzTKdLjk0R1VUaRg1l7cnlnTHvC/Cs+D/rFBvPE/GQW78tnRLsIL2co9ZFOgYiISBWLNx1gUHxYlbH+cSFs3Z9DcanDo1jTlidx1ZldKgtZgEB/Hy4e0oHJC3Z5FOtAViGFxeX0aBZUZfy0VqEs3JDiUay67kBmIZcP71xlbOTJ8ZQ7XOxIzvEo1qKNqSTEBFUWsgC+douEmCCPf29FDjfpRQ76tQiuMj64dRibM0o8igWwJ6eU09uGV5kdjgj0oWNUAIv253sUKzGnjM4tIyoL2V9ccWYXtuWUe5xbXeVwGXakF3HRkA5Vxsee3olNBwu9lJUcj82HSji1VWhlIQsQ5GunX4tgth7y/HVVV23LLOGkZkGVhSyAZVkMig9jcwN6nHJiaWZWRE4YtRXXD03CAjhU5CAy8NePiOwSJ/6+Nvx9PWvVbBYRxK6U3CPGd6Xm0i7Ws7bg8GB/ypxuihzuKi2jh4ocNI0I9ChWXefv58Peg/m0j/t1piIzrwSn203zyOA/ueeRoiMCySp1HzGeVeaiaaRnvzd/u4XNgtxSV5XnR3qRgwh/z9t4A+w20gurniAxxnCoyMHQNp49P8L87aRkFuJyubH/pnDfm5ZPmF/DOXdvt0GQv52UjELaxPx60mnvwXwig3y9mJlUV4S/nUNFR54oPFTkYEDLhtMeH/4njzMioOE8TjmxGs67u4iI1IjbLurNxz9nkl/mBKDE4eb9jRncMLp7lSLhWDx700C+X5rIjyv3AhWFylfzd7J8cxpP3zDQo1ihQX5cMKgdH2zIpMxZUZxlFTv4bEs2t1/Ux6NYdd2o/q25+42FpGcXAVBc6uC2V3+iS6tIosI8uyby0qEd2ZpZwvLkAowxGGNYmVLApvRixh5uUz1WvnYbw9qEMWFtOqWHj0FuqZOP1h9iVAfPWwQv7taEKduz2ZNdClRcO/jttmxKnG7ObBf2F/euKj7cn6aBdsZ/uBynqyK3rXuzeOHTVZzRpuEsEGazLM5sF84dr86joLhixjkrr4S7X5vPmW08+51J3XB623AW7Stg86FioOJ9cm5iLumFDvrGhng5u5rTo1kQJU7D9J05uA+v2bMjs4Q5e3IZ0a7hrHkgJ5YWgBKRE0Yzs/WD22148O1FfDBjMy0jAkjJLeO809oy4cEz8a/GCsRvfLuBx/+3jIhQf1wuQ1Gpg9fvGlatlXkLS8q57rmZzF+fQmyYP6l5ZdxzaR+euG5Ag1r10+12c8Z937Jy20HaxYazLz2fFtEhLHnjUqIjgv46wO+s3HqQq5+egdPhxAIsHzufPD6aU3u08CjOtAmzKHO6eXv1QdYeKCQm1I8D+eWM7hjJ1T2jq3UM3lubzpw9uUQF+lBY7sJmWTw8qAXdm3k2Aw0VHQSvrU7nQKGD5hGBpGQUcvVJ0ZzRwL4oO1yG9zYcYkVKIe1jw9h9IJ/hbcO45qRo7LaG8zpoTNalFfLmqoME+9oocxoCfGw8cFoL4sMbxvWyvzhQUM4ryw6QX+ok2M9OXpmLW/o1Z2CrhnPCSWrGsS4ApWJWRGqNitf6LSuvhJ0pubSJCSO2ieeFxW+VlzuZvGg3vnYbFw/t4NECRkeTfKiAlIxCurWOIjykYX3Z+62ktDzmrtlP307N6Ne5+XHFcrsNmxIzMcbQq31TbNUoeqZNmFX596xiBxnFTuJC/QitRovxbxWXO1maXEhkoJ2EFsf/pfZAQTkFZS7aRPjj34C3/MgpcZJe5KBFqC9h/rpyrL5zuQ2JOaX42W3Eh/s1qBN0v2WMITm/nFKnm3aRAfjoBIwchYpZEfE6FbMiDctvi1kREZHacqzFbMM9XSkiIiIiIiINlnpSRKRGaTZWRERERE4EzcyKiIiIiIhIvaNiVkREREREROodtRmLyHFRW7GIiIiIeIOKWRFp0VfrrQAAIABJREFU8LLySpi7Nhl/PzsjE+IJCvD1dkr1wqY9GWxKzKJDXDindI2pM9tEOJ1OrnhmJtv2ZXPpsI6Mv26gt1OqVO5wMXftfnILyxjWuyUtokO8nVKtMMawPbOE9EIHbSMDaB1xfNsjZZc42ZxeTJCvjV4xwfja68ZzTURE6jYVsyLSoH0wYwsPvL2IIb1aUlzm4JZX5vHZ42cxol+8t1Ors0rLnVz5zEzW7khn0EktWLvzEDFRQXz37LlEhgZ4NbcvftrJTf+cSai/nbYRAbzy2Rpe/WIt+ybfRESId3PbsDuD8x/5gdbNQ2keFczfX1vAA2P78o+r+3s1r5qWX+bk2UWpFJa5aB8VwMcbM+gaHci9A1tUqwj9dmsWk7dm0bN5EHllLt5efZBHhrSkQ5R3j6eIiNR9KmZFpMHavj+bR95byop3LqdTq0gAFm5I4dLx09nz+fWEBvl5OcO66flPVwOw+7Pr8PWx43Yb7nxtPve/vYgPHh7p1dxueXE2oztGck2vpliWRZnTzePz99Pj2k9I+eZmr+Xlcrm5dPx0Xrp1MGOHdwLgYHYRg+78igHdYzm9Tyuv5VaTxowbxeXjp9GlaSDXHz4GDpebfy5J5bvtWVzWPdqjeFsOFfPj7hzeOLstTYIqOiaWJefzwuIU3j23PXabZmhFROSPaQEoEWmwPp+7g2tHdassZAGG9m7JwB6xTF2W6MXM6rZPZ2/j6RsG4OtjB8Bms3jmhoF8vWAXDqfLa3kt2pBCmdPN5T2iK1ue/X1sXNOrGfkFpV7LC2DZljTCgvwqC1mAmKhg7rmkD5/M3ubFzGpWSZmT6Sv2MrZbk8pj4Gu3cUWPaOYn5Xscb/7ePMZ0iqosZAFObRVGeIAPWzKKayxvERFpmFTMikiDVVTqICLkyNnXiGB/ikqdXsiofigucxIeXPUayJBAX5wug8ttvJQVHMguxm7jiFbWYF8bxntpAb881468bjQ8xJ+ikobzXHM4XRgDAT5Vvz6E+Nkpc7k9jlfmNAT7HvlVJNjXRpnTywdVRETqPBWzIuKRsXmvVPlTl509oC2fzN5Ocamjciwtq4gZK/dyVv/WXsysbht9Shvem7a5ythHM7cxqEcsAX7euzrlkiHtsFkWy5MLqozP2JWDz1EKohNp0Ekt2JSYyda9WZVjTpebD6Zv4ewBbbyXWA0LC/anZ7smLNibV2V81u5cEmI9X+yqX4tg5ibmVTlJkpJfxu7sUro3CzzufEVEpGHTNbMi0mCd3qclA7rHMPD2L7l+dDeKy5xMmPozD17ej1bNQr2dXp319A0DGXb3ZHan5nJ6n1as2Z7O1GWJzHz5Aq/m5ePjw3lDOvCfBbvYmF5M+8gAlifnsy2jhK+eHuPV3EIC/Xj1jiGMuO8bbh5zEjGRQUyau53wEH+uPKOzV3OraW/eN4KzHviW3bnltAnzY8PBIhJzSvnnGZ6fIBocH8biffk8PHcfp7cJJ6/MyazdudzQtxlBvvZayF5ERBoSy3i7N8tD7br1NM9NmuHtNEQarbo+G/t7xhh+XLmXqcsS8ff14fIRnRjQLdbbadV5uYVlTJy5lU2JmXSIi+D60d1oHhXs7bQA+HjmFh54ezFOh5OI8CC+f+5cerZv6u20ANiclMnEWdvIKyznjIR4LhzcHh97w2uCSssq4oPpm1m0cBttI/0Z3ja82sWny21YmVrA+rQignztnN4mjDaRWslYRKQxO//z7WuNMQl/dTsVsyLyp+pb8SoiJ860CbO8nYKIiDRAx1rMNrzTxSIiIiIiItLg6ZpZEalCM7EiIiIiUh9oZlZERERERETqHRWzIiIiIiIiUu+ozVikkVNbsYhU15hxo6r8rAWhRETkRNLMrIiIiIiIiNQ7mpkVEWkAiksdPPXRSibO3kphiYOz+rfh+ZtPpWPLSG+nRlpWEY+8t5QpS/Zgt1lcOqwjz910GlFhnu8laozhv99v4rVvNrA/vYCEzs148vqBDO/bqhYy967Zq/fx+HtL2ZSURavoYB644mRuPrcHlmV5OzURkWopc7r57OdM5iflUeZy0zc2hGt6NSU21M/bqUk9pWJWpJFRW3HDdOUzM/HztbH4jcuICg3g/embGX7vN6x770qaRgR5La/Scidn3PcN557ajm0Tr8HpcvP8p6sZ/dB3LHtrLHa7Zw1CL3+xli9/2slH/xjJSW2j+XHVXq569kcmPzWG005qUUuP4sSbvz6Zq5+ewc29m/JQ7w4k5pTy8sTlFJc5uOfSvt5OT0SkWl5emoqfj40Xz2xNiJ+d2XtyeXTefv4zui1h/nZvpyf1kNqMRUTquc1JmazbeYhJj51Fh7gIosICePCKBEad3JoPf9zq1dy+Wbibls1CeeGWQcREBdOyaShv3XM6xsDsNfs9ilXucPGvL9fx5ZNnM6BbLMGBvlwytCPP3ngqL3+xtpYegXc8//EK/tajCQNbheLvY6Nr0yDuPSWGFz5djdPl9nZ6IiIeS8opJSm3jPsGtiA21I9QfzsXd2tCr5hg5ibmejs9qadUzIqI1HNbkrIZ0D0GX5+qZ7UH94xjS1KWl7KqsHVvFkN6xlUZsyyLwT3j2LrXs9zSc4rx87XRIS6iynh1YtV12/bn0L1Z1Rn1+HB/HE4X2fmlXspKRKT6kvPK6RwdiI+t6qUS3ZsGkpxX5qWspL5Tm7FIA6e24oavS3wkq7en43S58flN2+7yLWl0iffuNbOd46P4cv6OKmPGGJZvSeORv53sUaxmEYGUlbtISsujbWx45fjyLWl09vLjrGmdWkawPaOEZsG+lWMp+WXY7TYiQ/29mJmISPXEhfmxa1MJLrfB/puCdntWCXFhel+T6tHMrIhIPderQ1O6to7i+hdmk5JRQHGpgze+2cDUZYnccHZ3r+Z26bCO7ErO5ckPl5NbWEZmXgn3v72I0nIno/u38SiWv58Pd13cmyue/pENuzNwudxMW57II+8t5YGx/WrnAXjJP645hY9/zmTtgULcxrAnu5T/rErngcv7HTEDLyJSH7SPCqBFqD+vr0wjq9hBmdPN1B3ZrEkt5Ix24X8dQOQoLGOMt3PwSLtuPc1zk2Z4Ow2RekMzs41DYUk5j763jImzt1Fc6mTkyfG8eMsgurVp4u3USD5UwIPvLOaHpYnYDq9m/NKtg6q1MJXbbXht8npe+2Y9KRmF9O3YjCevH8DZA9rWQube9cPSRJ54bwmb9+cQFxXEfZcncNfFvev0asbaZ1ZE/kyxw8WnGzOYvzefMqebvrHBXNu7Ga3CNTMrVZ3/+fa1xpiEv7qdilmRBk7FbONijMEYsNnqXsHzy+dNTRVjbrepk4+zptWnx6liVkSOhTEGA9jq8Mk58a5jLWZ1zayISANiWRZ19btBTc8o1pcC73g1lscpIo2HZVnonU1qgq6ZFRERERERkXpHM7MiDZBai0VERESkodPMrIiIiIiIiNQ7KmZFRERERESk3lGbsUgDoLZiEakLxowbVeVnrW4sIiK1qdaKWcuyWgETgeaAASYYY1773W2GAd8DSYeHvjXGPF1bOYnIr8odLiYv3MWSTanENQ3h2rO60bJpqLfTAiCvsIxJc7azdV82XeIjuXpkVyJCtAfdX9m6N4vP5u6goLics05pw6iTW1d7JdzCknImzdnBz4mZdIgL528ju9IkPLBasYwxLPn5AN8u3IXdZjF2RGdO7hJTrVhSN+xIzuGhdxaxLy2f3h2b8dJtg2kW6fm+wY3N7uxSluzLx2kMA1uG0q1pYJ3eN1hEpK6rzTZjJ3C/MaYbMAC4w7Ksbke53WJjTO/Df1TIipwARSUOhtz5JS9/uAxX6iFWLttJ7+s/Zf76ZG+nxr6D+fS5aRILN6bSJT6SZVvS6H3jp+xJzfV2anXahz9uYfi93+A2hlbNQvnHhCVc/dxM3G7P9xI/kFlIv5s/Y9bqfXSJj2T9rgx63TiJbfuyq5XbA28t5Krx08jancrBHSmc99AUnpu4slqxxPsmL9hJnxs+JT8tmwFNfPh5czIdLv+ArUlZ3k6tTvt2axbPLUrB38ci3N/O6yvT+N+6Q95OS0SkXrN+2cS+1v8hy/oeeNMYM+c3Y8OAB4wxY441TrtuPc1zk2bUQoYi9cfxthX/89NVzJy3lQcGxFTOCqw9UMjEbbns+vwGr+5recXTP9KtdRSPX3tK5diLn61m1fZ0vnn6mN8qGpXcwjLaX/Ehy94aS+f4SABKy50MuvMrHrvmFC4Y1N6jeDe9NIfo8EBeuGVQ5dib325g2vIkZr58oUex1uxI57yHvuPVM+MJ8bMDkFPi5N7Z+1gx4Uo6xEV4FE+8L+a8/3JJl0jO6vDrsZuwNp395RYbPvxblduqzbhCemE598/ax+tntyUqsKIprqjcxT0zk3jg1Dg6R1ev60FEpKE6//Pta40xCX91uxOyAJRlWW2APsDRTsUPtCxro2VZP1qW1f1E5CPS2E1ZtItRbcOqtLf1jQ2mvNzJ9v3Vm32rKdOXJ3HbBT2rjN1+QS+mL0/iRJ18q29+WpfMgO4xlYUsQICfDzee3Z2pSxM9jjd9RRK3nV/1GNw0pgcLN6ZSWu70KNbUpYkMahlSWcgCRAb6MKBlKDOWJ/3JPaUuSskoIKeojDPahVcZH9Mpkt3qnvhDa9OK6B8XUlnIAgT72RnWJpxVqYVezExEpH6r9QWgLMsKAb4B7jHG5P/uf68DWhtjCi3LOhuYAnQ8SoxxwDiA6Ji4Ws5YpO6p6QWeAvx8KHG6q4y5DZQ73QT4eXdduAA/O0UlDqJ/c31mYYkDf1/7n9yrcQvws1NccmSRWVDiINDf8+MZ4OdDUamjylhJmRO7zcLu4ax9gL+dMteRJyFKXQZ/Lz/XxHMBvnYMUO5y42P79TVZ4nB7/NxoTPxsFqUu9xHjJU43oX56bxMRqa5anZm1LMuXikJ2kjHm29//f2NMvjGm8PDfZwC+lmVFH+V2E4wxCcaYhNDIqNpMWaRRuGZ0d6bszKX0NwXtzN25tIkNp12L8D+5Z+27YkRnxn+4ovJaT2MM4z9czpVndNZCKX9gRN9W7EjOYfbqfZVjB7OLeOf7TVwxorPH8a4Y0ZmnPlqJ8/CXb2MMT3+8kgsHt8fXx7Mv3mNP78yS5AJS8ssqxxKzS1mXVshFQzxrfxbvi44IomWTYL7cnFXZKeF0Gyb9nEm/zs29nF3d1b9lKBsPFrErq6RyLK2gnIV78xncOsyLmYmI1G+1uZqxBbwPbDPGvPoHt4kB0o0xxrKs/lQU11pBQqSWXT+6O8t+TuWOHxPpExtMWqGDQhfMevVib6fGszedyoWPTaX7dRM5rUcLlm9Jo2lEIFOeO8/bqdVZ/n4+fPnk2Vz25HR6d2hGVKg/s1bv48HL+3HaSS08jvfY3/pz6ZPT6XrNRIb0jGPtznQC/XyY+sL5Hsdq1yKcV+8cxj1vLKBXbDAut2HLoWLef3gkTSO0+m199MML5zPkzq9YlVpIx+hANqQVEhLsz6bnzvV2anVWmL+du06J5ckFyXRrGoSPzWLjwSKu79OMFqF+3k5PRKTeqrUFoCzLGgQsBn4Gfpn+eQSIBzDG/NeyrDuB26hY+bgEuM8Ys+zP4moBKGkMTtS+sVv3ZrFsSxqxUcGM6t8aH/sJuYz+LxljWL4ljS17s+naOpLTerTQrOwxKC51MH3FXgpLyjkzIf64tloyxrB6ezob92TSsWUEQ3vFHdcxyMwrYcaKvdhtFucMbKutluo5p9PNa9+s5+fETEb0i+dvI7se9XZaAKqqonIXaw4U4jLQLzaY8AC12ouIHM2xLgB1wlYzrikqZqUxOFHFrIhIbVIxKyIi1VGnVjMWERERERERqUnqbxGpAzQTKyIiIiLiGc3MioiIiIiISL2jYlZERERERETqHbUZi4iISK0YM25UlZ+1IJSIiNQkzcyKNGIlZU427M7gQGaht1ORGrJgQwqTF+yktNzp7VREREREapVmZkUaqXe+38T4D5YT2ySYA1lFDO0Vx/sPnUm49v+sl5b8nMpl46dTWu4iJNCXG1+ayz2X9uGp6wd6OzURERGRWqFiVsQLvL168cxVe/nXl2tZ9MaldImPorjUwX1vLeKml+fy9VPneDU38ZzT6eb8R6by0BX9uH9sP3zsNhZtTGXMP77n1O4tGNW/tbdTFBEREalxajMWaYTe+X4TT1w7gC7xUQAEBfjy6h1DmL8+mYPZRV7OTjz1zg+baBIWwENXJOBjr3hbH9IrjtvOO4lnPl7h5exEREREaoeKWZFGKD27mA5x4VXGggJ8iYkK5lBOsZeykupKPJBHh7hwLMuqMt45PorcwjIvZSUiIiJSu9RmLHICeLut+PdO69GCbxft5tQeLSrHtiRlkZVfQudWkV7MTKrj8uGdOOP+b8kpKCUyNAAAYwyfzd1O387NvZydiIiISO1QMSvSCN13WV9OveNL3MZw4eAO7ErJ5dmJK3n6hlPx99PbQn1zSrdYTmrbhNPu+JKnbjiVJmEB/Pf7TWzcncEXT4z2dnoiIiIitULfWkUaobimISx7ayz//nod//fuEmKaBPPuAyM4M0ELBdVXS968jLteX8CD7yzC4XTTo10069+/iuiIIG+nJiIiIlIrVMyK1IK61lZ8NHFNQ3jl9iHeTkNqiM1m4817hvPmPd7OREREROTE0AJQIiIiIiIiUu+omBUREREREZF6R23GIjWkPrQWi4h405hxoyr/Pm3CLC9mIiIiDYFmZkVERERERKTe0cysSDVpJlZERERExHs0MysiIiIiIiL1jmZmpd5zOZ389O0kls/6AafTScKwkYy6/Ab8AwO9nVqj4XYbPvhxC5/N2U5puYtzBrblrot7Exrk5+3UGg1jDBNnbeOT2dsoLHFwVv/W3HNpXyJC/L2dmjRAS38+wGvfrCfpQD69OkTzwOX96BIf5e20RESkkVExK/XeO0/cQ05GOufdcCe+fv7M/uJDXrzzah5990vsPjX3FFdb8R+79dV5bN2bzf9dlUBooB///eFnzrz/Wxa+dgn+fnqbORHue2sRyzYf4JGr+xMVFsD707dw+t2TWfLmZQQH+no7PWlAvl+yh9v//RNPXDuAvpc3Ze6a/Qy7ezKzX7mQnu2bejs9ERFpRPQtU+q1xK2b2LlhDa98twA//wAAuvYbwJPXX8i6RXM4efhoL2fY8G3fn820ZUnsmnRdZdE0pFccI+//lq8X7ubqM7t4OcOGb9/BfD6ds53dk64j/PBM7KCTWnD+o1P5dM52bjnvJC9nKA2FMYZH3lvKxEdGMaJfPAAnd4khJNCPZyau4uunzvFyhiIi0pjomlmp13ZtWkevQadXFrIANpuNhGEj2blxrRczazyWb0njzIT4KrN/lmVx/qD2LN98wIuZNR4rtx1kaK+4ykIWKo7BBYPas0zHQGpQbmEZqZlFDO/bqsr4BYP1XBMRkRNPxazUa5FNm5G2N/GI8QN79xDRtJkXMmp8YpsEszM554jxHck5xDYJ9kJGjU9sk2B2peRijKkyvjM5hxbROgZSc4IDfLHbLFIzC6uM70jOoUWTEC9lJSIijZWKWanX+gwewaHU/cz5eiJulwtjDKvm/ciGJT8x6OyLvJ1eo3BGv3hyi8r491frcLrcGGOYsSKJr+bv4tqzunk7vUbhtB4tsNssnv90NeWOitfBvLX7+fDHrdx4Tg9vpycNiJ+vnetHd+O2V38iO78UgL0H83nw7cXcfmFPL2cnIiKNjfX7M/l1XbtuPc1zk2Z4Ow2pQ1KTdvPu+PvISEvBx9eXgMBgbn7iJTr1SqjRf0cLQP2xxAN5XP/CbHYk5xDk74Ofr5137hvO6X1a/fWdpUYkHyrghhdmsykxi9CgipbvN+85nbP6t/FuYtLglDtc3PvmQj6ft4PYJsEcyi3h/sv68vCVCViWdcxxpk2YVYtZiohIfXb+59vXGmP+8su8illpMDIOJON0OIiJb+vRF6o/ouLVc/vT8ykpc9GpVUSNHAPxXPKhAopKHXRqGYnNpmMgtSenoJQDmUW0iQmr1orZKmZFROSPHGsxq9WMpcFo2kKzgN4W3zzM2yk0eq2ahXo7BWkkIkMDiAwN+OsbioiI1BJdMysiIiIiIiL1jmZmRQ5TW7GIyIkzZtyoKj+r7VhERDylmVkRERERERGpd1TMioiIiIiISL2jNmNptNRWLCIiIiJSf6mYFRERqSGZucU8M3EVWfmlXD+6GyP6xXs7pUqZeSVMXZaIMTBmYFuaRQZ5OyUREZHjojZjERGRGvDe1J+Jv+R/LFy2i7TEg5z/f98z6PYvcbvd3k6NSXO20/nqj5m9ej9z1+6n6zUT+XjmVm+nJSIiclw0MyuNhtqKRaS2lJY7uef1BTx4Whwnx4UAUNCnGffMTOL5Sat57G+neC23lIwC7n5jAYvfuJRubZoAsGN/DoP+/hWn92mp/aFFRKTe0sysiIjIcfrvDz/TNNi3spAFCPW3c0m3Jnzi5RnQbxfu5oJB7SsLWYDO8ZFcMrQDkxfu9mJmIiIix0fFrIiIyHEqKXPiZ7eOGPezW7hcxgsZ/arc6SbQ/8hGrEB/H8odLi9kJCIiUjNUzEqDNTbvlSp/RERqy83ndCc5r4w92aWVYw6X4YedOYwe2NaLmcGYU9syeeEu0rOLKscycov58qednHtqOy9mJiIicnx0zayIiMhxio4I4s6LevOP7zZwRrsIogJ9mJuYh2+AL/+6bbBXc+sSH8WdF/bm5Fu+4LrR3bCAj2Zu5eYxPejetslf3l9ERKSuUjErDYpmYEXEW16+fQjnntaOpz9eycGiMm6/tB8PXZGAj4/3m6Ae/Vt/Rp/ShskLd2GM4ZtnxpDQubm30xIRETkuKmZFRERqyJBeLZn7aktvp3FUfTs1o2+nZt5OQ0REpMZ4/3SxiIiIiIiIiIc0Myv1mtqKRUQahjHjRlX5edqEWV7KRERE6gvNzIqIiIiIiEi9o2JWRERERERE6h0VsyIiIiIiIlLvqJgVERERERGRekfFrIiIiIiIiNQ7Ws1Y6hWtXiwiIiIiIqCZWREREREREamHVMyKiIiIiIhIvaM2Y6nT1FYsIiIiIiJHo5lZERERERERqXdUzIqIiIiIiEi9ozZjqVPUViwiIiIiIsdCM7MiIiIiIiJS76iYFRERERERkXpHbcbiVWorFhGRoxkzblSVn6dNmOWlTEREpK7SzKyIiIiIiIjUOypmRUREREREpN5Rm7GcUGorFhERERGRmlBrxaxlWa2AiUBzwAATjDGv/e42FvAacDZQDFxnjFlXWzmJSO1ZtvkAn83dQWm5k3MGtuW8U9tht6v5Q0RERERqR23OzDqB+40x6yzLCgXWWpY1xxiz9Te3GQ10PPznFOCdw/+VBkIzsY3Dy5+v4a0pG7n9gl6EBvry/Cer+fKnnXz2+GhsNsvb6YmIiIhIA1RrxawxJg1IO/z3AsuytgFxwG+L2fOBicYYA6ywLCvCsqzYw/cVkXogNaOQFz9fw6YPrqZFdAgAN5zdnQG3fcnMVXs5e0BbL2coIiIiIg3RCekBtCyrDdAHWPm7/xUHJP/m55TDYyJST8xbl8yZCa0rC1kAfz8f/jaqKz+u3Ou9xERERESkQav1YtayrBDgG+AeY0x+NWOMsyxrjWVZawpysms2Qfn/9u491rK6ugP4d4Xh/dRKWi0q1CoN0AZwAgXRWqUiOAGNttXaP3wko401aY1t+tZGrPFRWyOJZESqjRRTCTRmtIxWimI1UpDXIIKhUgrSQgsiWh5FVv84m/YOMDAz9567Z9/7+SSTuXufvX/3eyB7btZda+8Di7LPnrvmznvue9T+O79/X/bda7cREgEAsBrMtZitql0zK2TP6e7zH+OQW5M8fcH2QcO+LXT3hu5e291r933Sk+cTFtghJx97cK6+8Y5ceOlN/7fv27fclY99bnNee+LPjBcMAIAVbZ5PM64kH0tyXXd/cCuHfSbJb1bVpzJ78NPd7peFadlz9zU570/X5Vfe+dk8+6ADss+eu+Vr196W9735hBx+yI+NHQ8AgBWqZs9emsPCVSckuSTJNUkeGnb/QZJnJEl3nzkUvGckeWlmH83z+u6+7PHW/anDfq7ffc7n5pKZpeEJxqvT/Q88mIuuuCX33v9gXnT003PAPruPHQlYwTZu2DR2BADm5LRzv3V5d699ouPm+TTjryR53M/kGJ5i/JZ5ZQCWz+67rcnJxx48dgwAAFaJZXmaMQAAACyluXVmWT2MFQMAAMtNZxYAAIDJUcwCAAAwOcaM2W7GigEAgLHpzAIAADA5ilkAAAAmx5gxT8hYMQAAsLPRmQUAAGByFLMAAABMjjFjHsVYMQAAsLPTmQUAAGByFLMAAABMjjFjjBUDMDnr1p+0xfbGDZtGSgLAWHRmAQAAmBzFLAAAAJOjmAUAAGByFLMAAABMjgdArUIe+AQAAEydziwAAACTo5gFAABgcowZrwLGigEAgJVGZxYAAIDJUcwCAAAwOcaMVyBjxQAAwEqnMwsAAMDkKGYBAACYHGPGK4TRYgBWs3XrT9pie+OGTSMlAWC56MwCAAAwOYpZAAAAJseY8UQZKwYAAFYznVkAAAAmRzELAADA5BgznghjxQAAAP9PZxYAAIDJUcwCAAAwOYpZAAAAJkcxCwAAwOQoZgEAAJgcTzPeSXl6MQDsuHXrT9pie+OGTSMlAWBedGYBAACYHJ3ZnYROLAAAwLbTmQUAAGByFLMAAABMjjHjkRgrBgAA2HE6swAAAEyOYhYAAIDJMWa8TIwVAwAALB2dWQAAACZHMQsAAMDkGDOeE2PFAAAA86MzCwAAwOQoZgEAAJgcY8ZLyGgxAOyc1q0/aYvtjRs2jZQEgKWiMwsAAMDkKGYBAACYHGMF2OPkAAAIIElEQVTGi2CsGAAAYBw6swAAAEyOYhYAAIDJUcwCAAAwOYpZAAAAJkcxCwAAwOR4mvF28PRiAACAncPcOrNVdXZV3V5Vm7fy+gur6u6qunL48yfzygIAAMDKMs/O7MeTnJHkrx/nmEu6e90cMyyKTiwAAMDOaW6d2e7+cpI757U+AAAAq9fYD4A6rqquqqq/r6rDR84CAADARFR3z2/xqoOTbOzuIx7jtf2SPNTdP6iqU5J8qLufvZV11idZP2wemuT6+SRetKck+c+xQ8DIXAfgOoDEdQCJ62BHPbO7D3yig0YrZh/j2JuSrO3uyf7PrqrLunvt2DlgTK4DcB1A4jqAxHUwb6ONGVfVT1RVDV8fM2T5r7HyAAAAMB1ze5pxVZ2b5IVJnlJVtyR5R5Jdk6S7z0zyqiS/UVUPJrk3yat7nm1iAAAAVoy5FbPd/ZoneP2MzD66ZyXZMHYA2Am4DsB1AInrABLXwVzN9Z5ZAAAAmIexP5oHAAAAtptidolU1S5VdUVVbRw7C4ylqm6qqmuq6sqqumzsPDCGqjqgqs6rqm9V1XVVddzYmWA5VdWhw8+Bh/98v6p+a+xcsNyq6rer6tqq2lxV51bVHmNnWmmMGS+RqnpbkrVJ9uvudWPngTGshI/YgsWqqk8kuaS7z6qq3ZLs1d3fGzsXjKGqdklya5Jju/tfx84Dy6WqfjLJV5Ic1t33VtXfJvlcd3983GQri87sEqiqg5K8LMlZY2cBYDxVtX+SFyT5WJJ09wMKWVa5Fye5USHLKrUmyZ5VtSbJXkm+O3KeFUcxuzT+MsnvJnlo7CAwsk7y+aq6vKrWjx0GRnBIkjuS/NVw68lZVbX32KFgRK9Ocu7YIWC5dfetST6Q5OYktyW5u7s/P26qlUcxu0hVtS7J7d19+dhZYCdwQncfneTkJG+pqheMHQiW2ZokRyf5SHcfleSHSX5v3EgwjmHM/tQknx47Cyy3qnpSktMy+yXn05LsXVW/Pm6qlUcxu3jPS3LqcK/gp5K8qKo+OW4kGMfwW8h09+1JLkhyzLiJYNndkuSW7v76sH1eZsUtrEYnJ/lGd//H2EFgBCcm+U5339Hd/5Pk/CTHj5xpxVHMLlJ3/353H9TdB2c2SnNRd/utC6tOVe1dVfs+/HWSlyTZPG4qWF7d/e9J/q2qDh12vTjJN0eMBGN6TYwYs3rdnOTnq2qvqqrMfh5cN3KmFWfN2AGAFePHk1ww+/c6a5L8TXdfOG4kGMVbk5wzjFj+S5LXj5wHlt3wS81fSvKmsbPAGLr761V1XpJvJHkwyRVJNoybauXx0TwAAABMjjFjAAAAJkcxCwAAwOQoZgEAAJgcxSwAAACTo5gFAABgchSzALDEquqFVbVxW/cvwfd7eVUdtmD74qpauw3nPXUp8lTVgVXlo7gAWFaKWQCYvpcnOewJj3q0tyX56GK/eXffkeS2qnreYtcCgG2lmAVg1amqvavqs1V1VVVtrqpfHfY/t6q+VFWXV9WmqnrqsP/iqvpQVV05HH/MsP+YqvpaVV1RVV+tqkO3M8PZVXXpcP5pw/7XVdX5VXVhVX27qt634Jw3VtUNwzkfraozqur4JKcmef+Q71nD4b88HHdDVT1/KzFemeTCYe1dquoDw/u7uqreOuy/qareM6x9WVUdPfy3ubGq3rxgrb9L8tptff8AsFhrxg4AACN4aZLvdvfLkqSq9q+qXZN8OMlp3X3HUOC+O8kbhnP26u4jq+oFSc5OckSSbyV5fnc/WFUnJvmzzArEbfGHSS7q7jdU1QFJLq2qfxheOzLJUUnuT3J9VX04yY+S/HGSo5Pck+SiJFd191er6jNJNnb3ecP7SZI13X1MVZ2S5B1JTlz4zavqkCR3dff9w671SQ5OcuTwfp684PCbh/f+F0k+nuR5SfZIsjnJmcMxlyU5fRvfOwAsmmIWgNXomiR/XlXvzawIvKSqjsisQP3CUAzukuS2BeecmyTd/eWq2m8oQPdN8omqenaSTrLrdmR4SZJTq+rtw/YeSZ4xfP3F7r47Sarqm0memeQpSb7U3XcO+z+d5DmPs/75w9+XZ1akPtJTk9yxYPvEJGd294PD+7xzwWufGf6+Jsk+3X1Pknuq6v6qOqC7v5fk9iRPe/y3DABLRzELwKrT3TdU1dFJTklyelV9MckFSa7t7uO2dtpjbL8ryT929yuq6uAkF29HjEryyu6+foudVcdm1pF92I+yYz+vH15ja+ffm1kBvT1rPfSIbA8tWHuPYU0AWBbumQVg1amqpyX57+7+ZJL3Zza6e32SA6vquOGYXavq8AWnPXxf7QlJ7h46p/snuXV4/XXbGWNTkrfW0AauqqOe4Ph/TvILVfWkqlqTLceZ78msS7w9bsiWHdsvJHnTsHYeMWa8LZ6T2dgxACwLxSwAq9HPZnaP6pWZ3U96enc/kORVSd5bVVcluTLJ8QvOua+qrsjsHtE3Dvvel+Q9w/7t7Z6+K7Ox5Kur6tphe6u6+9bM7sm9NMk/Jbkpyd3Dy59K8jvDg6Se9dgrPGq9Hya5sap+eth1VpKbhzxXJfm17Xs7+cUkn93OcwBgh1X3I6emAICFquriJG/v7stGzrFPd/9g6J5ekOTs7r5gEeu9Islzu/uPliDblzN7eNZdi10LALaFziwATMc7h27y5iTfyezjcHbYUAjftNhQVXVgkg8qZAFYTjqzAAAATI7OLAAAAJOjmAUAAGByFLMAAABMjmIWAACAyVHMAgAAMDmKWQAAACbnfwEsvHIsp4qd7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure(figsize=(16, 9))\n",
    "\n",
    "plt.pcolormesh(xx, yy, pred, cmap=plt.cm.Paired);\n",
    "plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=plt.cm.Paired);\n",
    "plt.xlabel(iris.feature_names[0]);\n",
    "plt.ylabel(iris.feature_names[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>\n",
    "**References**\n",
    "\n",
    "- Agresti, A. (1996). An introduction to categorical data analysis (Vol. 135). Wiley.\n",
    "- Friedman, J., Hastie, T., & Tibshirani, R. (2001). The elements of statistical learning. Springer series in statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression Excercise - Space Shuttle Challenger Disaster\n",
    "\n",
    "> In January of 1986, the space shuttle Challenger exploded shortly after launch. Investigation of the incident suggests that the cause of the crash maybe related to the rubber O-ring seals used in the rocket booster. O-rings are more brittle at lower temperatures and the temperature at the time of the launch was 31 degrees Farenheit.\n",
    ">\n",
    "> Can one predict the failure of the O-rings in advance? There is data regarding the failure rate of O-rings for 23 previous shuttle missions. There is 6 O-rings used for each shuttle mission, and the number of O-rings damaged for each mission recorded.\n",
    ">\n",
    "\n",
    "### Question: Assuming a multinomial model like above, what is the probability of no O-rings fails at 31 degree Farenheit? $\\mathbb{P}(y=0 | \\ X = 31$)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Farenheit</th>\n",
       "      <th>Number of O-rings damaged</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Shuttle Mission Number</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>53</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>57</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>58</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>63</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>66</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>67</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>67</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>67</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>68</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>69</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>70</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>70</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>72</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>73</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>75</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>75</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>76</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>76</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>78</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>79</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>81</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        Farenheit  Number of O-rings damaged\n",
       "Shuttle Mission Number                                      \n",
       "1                              53                          5\n",
       "2                              57                          1\n",
       "3                              58                          1\n",
       "4                              63                          1\n",
       "5                              66                          0\n",
       "6                              67                          0\n",
       "7                              67                          0\n",
       "8                              67                          0\n",
       "9                              68                          0\n",
       "10                             69                          0\n",
       "11                             70                          1\n",
       "12                             70                          0\n",
       "13                             70                          1\n",
       "14                             70                          0\n",
       "15                             72                          0\n",
       "16                             73                          0\n",
       "17                             75                          0\n",
       "18                             75                          1\n",
       "19                             76                          0\n",
       "20                             76                          0\n",
       "21                             78                          0\n",
       "22                             79                          0\n",
       "23                             81                          0"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load and print the O-rings data to gain insight\n",
    "link = '/home/akhil/machine-learning-basics/data/Orings.csv'\n",
    "df = pd.read_csv(link, header=None, names=['Farenheit', 'Number of O-rings damaged'])\n",
    "df.index.names = ['Shuttle Mission Number']\n",
    "df.index +=1\n",
    "\n",
    "#df.loc[1, 'Number of O-rings damaged'] = 1\n",
    "# df.loc[1, 'Number of O-rings undamaged'] = 5\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPd57/HPg5BAYHbEohmD8A5mHWy874Z4wdvIafbbJm7c3sRpmt62jpP2OmnixmmbtkmzOE6I47ax01xLeIsXnMQLTmJjGLEY8A7YI7EIs2MWIT33j3MUBiyJkXRGMyN936/XvGbmN2fO75nzkuaZ8zvP+R1zd0RERPrlOwARESkMSggiIgIoIYiISEgJQUREACUEEREJKSGIiAighCAiIiElBBERAZQQREQk1D/fAXTG6NGjvaqqKt9hiIgUlWXLlm1194pjLVdUCaGqqoqlS5fmOwwRkaJiZhuyWU5DRiIiAighiIhISAlBREQAJQQREQkpIYiICJDnKiMzWw/sBpqBQ+5+Rj7jERHpywqh7PQSd9+a7yAi0dgI69dDVRVUHLPkV0SkoGjIKCr33w8TJ8LcucH9/ffnOyIRkU7Jd0JwYJGZLTOzm/McS9c1NsJNN8G+fbBzZ3B/001Bu4hIkch3Qjjf3RPAlcBnzezCoxcws5vNbKmZLW0s1C/Y9euhrOzIttLSoF1EpEjkNSG4e314vwVYCMxpY5m73f0Mdz+jolDH5auq4ODBI9uamoJ2EZEikbeEYGaDzWxI62NgHvByvuLplooKWLAAysth6NDgfsECHVgWkaKSzyqjscBCM2uN4z53fyKP8XTPRz4Cl1+uKiMRKVp5Swju/hYwI1/950RFhRKBiBStfB9UFhGRAqGEICIigBKCiIiElBBERARQQhARkZASgoiIAEoIIiISUkIQERFACUFEREJKCCIiAighiIhISAlBREQAJQQREQkpIYiICKCEICIiISUEEREBlBBERCSkhCAiIoASgoiIhJQQREQEUEIQEZGQEoKIiABKCCIiElJCEBERAPq394KZ/Qfg7b3u7n8RRQBmVgIsBerdfX4U6xQRkc7raA9hKbAMGAgkgNfD20ygLMIYPg+sjXB9IiLSBe3uIbj7vQBm9r+B8939UPj8LmBxFJ2bWRy4GrgD+Kso1ikiIl2TzTGEEcDQjOfHhW1R+Hfgb4GW9hYws5vNbKmZLW1sbIyoWxEROVo2CeFOoM7Mfmpm9wIp4B+727GZzQe2uPuyjpZz97vd/Qx3P6OioqK73YqISDvaHTJq5e73mNnjwFlh063uvimCvs8DrjWzqwiOUww1s/92949HsG4REemkY+4hmJkBlwMz3P0hoMzM5nS3Y3e/zd3j7l4FfBj4jZKBiEj+ZDNk9H3gHOAj4fPdwPdyFpGIiOTFMYeMgLPcPWFmdQDuvt3Moiw7xd2fAZ6Jcp0iItI52ewhNIUnjzmAmVXQQVWQiIgUp2wSwneAhcAYM7sDeJ4IqoxERKSwZFNl9DMzWwZcBhhwvbvrzGIRkV7mmAnBzEYCW4D7M9pK3b0pl4GJiEjPymbIKAU0Aq8RzGXUCKw3s5SZzc5lcCIi0nOySQhPAVe5+2h3HwVcCTwKfIagJFVERHqBbBLC2e7+ZOsTd18EnOPuLwADchaZiIj0qGzOQ9hoZrcCPw+ffwjYHJaiqvxURKSXyGYP4aNAHHgwvE0I20qAP8pdaCIi0pOyKTvdCnyunZffiDYcERHJl2zKTisIrllwOsGspAC4+6U5jEtERHpYNkNGPwNeASYBXwXWAy/lMCYREcmDbBLCKHdfADS5+7Pu/ilAewciIr1MNlVGrWckbzSzq4EGYGTuQhIRkXzIJiF83cyGAf8H+A+C6yt/IadRiYhIj8umyujR8OFO4JLchiMiIvmSTZXRJIKy06rM5d392tyFJSIiPS2bIaMHgQXAI+jMZBGRXiubhLDf3b+T80hERCSvskkI3zaz24FFwIHWRndP5SwqERHpcdkkhGnAJwjOPWgdMnJ0LoKISK+STUL4IHCCux/MdTAiIpI/2Zyp/DIwPNeBiIhIfmWzhzAceMXMXuLIYwjdKjs1s4HAcwQX2ekPPODut3dnnSIi0nXZJIRcfUkfAC519z1mVgo8b2aPh1diExGRHpbNmcrP5qJjd3dgT/i0NLx5LvoSEZFjO+YxBDM728xeMrM9ZnbQzJrNbFcUnZtZiZktB7YAT7n7i1GsV0REOi+bg8rfBT4CvA6UA38KfC+Kzt292d1nElyic46ZTT16GTO72cyWmtnSxsbGKLoVEZE2ZJMQcPc3gJLwC/we4Ioog3D3HcDTba3X3e929zPc/YyKiooouxURkQzZHFR+z8zKgOVm9k/ARrJMJB0JL83Z5O47zKwcmAt8s7vrFRGRrsnmi/0TQAlwC7AXOB6ojqDv8cDTZraS4JKcT2VMtS0iIj0smyqjDeHDfQTXVI6Eu68EZkW1PhER6Z52E4KZraKDMlB3n56TiEREJC862kOYH95/Nrz/r/D+4+h8ARGRXqfdhNA6VGRmc909c2jnVjNLAV/MdXAiItJzsjmobGZ2XsaTc7N8n4iIFJFsyk5vAn5iZsPC5zuAT+UuJBERyYdsqoyWATNaE4K778x5VCIi0uOy2UMAlAhERHo7HQsQERFACUFERELZTH/9QTMbEj7+OzOrNbNE7kMTEZGelM0ewt+7+24zOx+4HFgA/CC3YYmISE/LJiE0h/dXA3e7+y+BstyFJCIi+ZBNQqg3sx8CHwIeM7MBWb5PRESKSDZf7H8EPAl8ILyQzUjgb3IalYiI9LhsEsJA4BngXTMbCRwguLpZ0Vjxzg6+8dhaXt20O9+hiIgUrGxOTEsRXBRnO2DAcGCTmW0GPh2eyVzQRg8ZQEk/40/uWcLIwWUkE3Gum1nJ6OMG5Ds0EZGCYe4dz2RtZj8CHnD3J8Pn84AbgZ8A33b3s3IeZeiMM87wpUuXdvn9zS3OC2+9S00qzVNrNjOnaiTJRJzLJo9hYGlJhJGKiBQOM1vm7mccc7ksEsIqd592VNtKd59uZsvdfWY3Y81adxNCpr0HDvHEy5uorUuzumEXV04dz42zYyQmjMDMIulDRKQQZJsQshky2mhmtwI/D59/CNhsZiVASzdizKvBA/pTPTtO9ew4DTv2sbCunr99YCWHWpzkrDjJRIzjRw7Kd5giIj0mmz2E0cDtwPlh028Jrq28E5jg7m/kNMIMUe4htMXdWZneSU0qzaMrN3JSxXFUz45x5bTxDB1YmrN+RURyKbIho0KS64SQ6eChFp5+dQu1qTS/e/NdLj51DMlEjAtOGk3/Ep2GISLFI8pjCKcAfw1UkTHE5O6XdjPGTuvJhJBp+96DPLqygZpUPfU79nHdjEqqZ8eZPH5oj8ciItJZUSaEFcBdwDIOT2NBPspN85UQMr3ZuIfaVJqFqXqGDSqjOhHj2pmVjBkyMK9xiYi0J8qEsMzdZ0cWWTcUQkJo1dLivLDuXWpT9SxavYnZE0eQTMSZO2WsSlhFpKBEmRC+AmwBFhKcpQyAu2/rZoDHA/8JjAWcYOK8b3f0nkJKCJneO3iIJ1dvojZVz8r0Tq6cOo5kIs6ZVSphFZH8izIhrGuj2d39hK4GF653PDDe3VPh9RaWAde7+5r23lOoCSHTxp1BCWttqp6Dh1q4YVaM6kScCaNUwioi+VF0VUZm9hDwXXd/qr1liiEhtHJ3VtXvpDZVzyMrGjihYjDJRJyrp6uEVUR6VrcTgpld6u6/MbNkW6+7e203Y8zsqwp4Dpjq7rvaW66YEkKmpuYWnnm1kZplaX775lYuOqWC6kScC05WCauI5F4UCeGr7n67md3Txsvu7p/qbpBhP8cBzwJ3tJVkzOxm4GaACRMmzN6wYUMU3ebNjvcO8sjKjdQsS5Pevo/rZlZSnYgzpVIlrCKSG5EMGZlZP+BGd/9FlMFlrL8UeBR40t3/9VjLF+seQnveatzzh+MNQwb2pzqchXXMUJWwikh0ojyovDSbFXWWBeU39wLb3P0vs3lPb0sIrVpanBfXbaM2lebJ1ZuYNWEEyUSMD5w+TiWsItJtUSaEO4GtwP8Ae1vbIyg7PR9YDKzi8CR5X3L3x9p7T29NCJneO3iIRas3U5NKs+KdHVw5dTzJRIwzq0bSr59KWEWk8wq+7LQr+kJCyLRp534eWl5PTSrNvqZmbpgVJzkrRtXowfkOTUSKSNGVnWajryWEVu7O6oZd1KTSPLKigYmjBlMdlrAOK1cJq4h0LCcJwczudvebuxVZN/TVhJCpqbmF515rpCaVZvHrW7nw5AqSiRgXnlJBqUpYRaQNuUoIKXdPdCuyblBCONLO95p4dFUDNcvSvL1tH9fOqCSZiHF65VBNmSEifxDlFdMybeliPJIDwwaV8rGzJvKxsyaybuteFqbS/Pl/L2NwWX+qZ8e4bmaMsSphFZEsZbWHEJ48hrvvyXlEHdAewrG1tDgvrd9GTSrNEy9vYuaEEVQnYsybMo7yMpWwivRFUZ2Y9hngi8BgwIDdwDfd/ftRBdoZSgids+9gM4vWBLOw1r29nSvCWVjnqIRVpE/p9pCRmf0dcC5wsbu/FbadAHzbzEa6+9cji1ZyoryshOtmBkNHW3bt58Hl9dz+0Gr2HDhEMhEjmYgzSSWsIhLqaC6jV4EZ7r7/qPZyYIW7n9ID8R1Bewjd11rCWpuq5+EV9UwYOYjq2XHmT6tk2CCVsIr0RlEcVPajk0HYuM/MWtp6gxQ+M2NqbBhTY8O47arTeO61RmpT9dz52CtccMpokrPiXHSqSlhF+qKOEkK9mV3m7r/ObDSzS4GNuQ1LekJpST8umzyWyyaP/UMJ613PvskXa1dyzYxgFlaVsIr0HR0NGZ0OPAQ8T3A1M4AzgPOA69x9dY9EmEFDRj1j/da91NbVs7AuTXlpCdWJONfPUgmrSLGKqspoIPBR4PSwaQ3ws7aGknqCEkLPamlxlm7YTs2yNE+s3sT0+DCqE3E+cLpKWEWKieYykkjtb2pm0ZrN1KbSpDZs5wOnByWsZ01SCatIocvVmcrSRw0sLeHaGZVcO6OSLbv289DyBr76yGp27z/EDbNiJBMxTqg4Lt9hikg3aA9BumV1w04Wpup5cHkDx48sJ5mIc8308QwfVJbv0EQkFMU1lX/t7peZ2Tfd/dbII+wCJYTCdai5hcWvb6UmlebZ1xo5/6TRJBNxLlYJq0jeRTFkNN7MzgWuNbOfE0xd8QfunupmjNKL9C/pxyWnjeGS08awc18Tj63ayN3PvckXaw6XsE6NqYRVpJB1tIdwI3ATcD5w9M9yd/dLcxzb+2gPofhseHcvtal6FtbVM6B/P5KJONfPqmT8sPJ8hybSZ0R5Cc2/d/evRRZZNyghFC93Z9mG7dSk0jy2ahPTYsNIJmJcMXUcg8pU2yCSS5GWnZrZtcCF4dNn3P3RbsbXJUoIvcP+pmZ+tXYzNcvSLNuwnblTxlGdiHH2CaNUwiqSA5GVnZrZN4A5wM/Cps+b2bnu/qVuxih91MDSEuZPr2T+9Eq27N7Pw8sb+Nov17LzvYPckIhxw6w4J41RCatIT8tmyGglMNPdW8LnJUCdu0/vgfiOoD2E3m1Nwy4W1qV5cHkDlcPLqU7EuGZ6JSMGq4RVpDuiPjFtOLAtfDysy1GJdGBK5VCmVE7h1itOY/EbW6lN1fPPT7zKuSeNIpmIc8mpYyjrrxJWkVzJJiF8A6gzs6cJSk8vJLiKmkhO9C/pxyWnjuGSU8ewa38Tj6/ayILn13Fb7SqumT6eZCLO9PgwlbCKRCzbg8rjgTPDp0vcfVMknZv9BJgPbHH3qcdaXkNGndDYCOvXQ1UVVFR0fZlC0djIO6vfpHbXQGpf3U5pST+SiRjXz4xRObyXlrCuXQtLlsCcOTB5cr6jObZi+nvqY7IdMspq/9vdN7r7w+EtkmQQ+ilwRYTrE4D774eJE2Hu3OD+/vu7tkyhCGM9/vor+PyHz+WZ2CbuTE7jnW37uPLbi/nYj1+gZlmavQcO5TvS6HzuczBlCvzJnwT3n/tcviPqWDH9PUm78j6XkZlVAY9qDyEijY3BP+S+fYfbysthw4bDv9qyWaZQHCPW1hLWhal6lqzfxtzJY6meHefsE0ZRUqwlrGvXBkngaGvWFOaeQjH9PfVRvWa2UzO7GbgZYMKECXmOpgisXw9lZUf+c5aWBu2t/5zZLFMojhFrZglr4+4DPLyigX98bC3b9h4MZ2EtwhLWJUvaby/EhFBMf0/SoQ4TQlhiutrdT+uheN7H3e8G7oZgDyFfcRSNqio4ePDItqamoL0zyxSKTsRaMWQAN50/iZvOn8Qrm3axMFXPR3/0AuOHDQxmYZ1RychiKGGdM6dz7flWTH9P0qEOjyG4ezPwqpnpp3mxqKiABQuCXfahQ4P7BQuO/KWWzTKFoouxnjZuKLddNZnf33YZfzXvVFJvb+eif36am/9zKU+8vIkDh5p76AN0weTJcMstR7bdckth7h1Acf09SYeyOTHtOWAWsATY29ru7tdGEoCOIeRGL6wy6m6su/c38fiqTdSk0ry2eTfzp1eSTMSYefzwwixhVZWRRCTKye0uaqvd3Z/tYmyZ674fuBgYDWwGbnf3Be0tr4QgUXln23s8WFdPbV09ZlCdiHP9rBix3lrCKn1a1JPbTQROdvdfmdkgoMTdd0cQZ6coIUjU3J3U2zuoTaX55aqNTBk/lGQizpVTxzF4QMHXXIhkJco9hE8TVPmMdPcTzexk4C53vyyaULOnhCC5tL+pmd+8soXaVJoX1wUlrMlEnHNOLOISVhGiLTv9LMFspy8CuPvrZjamm/GJFJyBpSVcNW08V00bz9Y9B3h4eQPffOIVGncf4PpZMaoTMU4eOyTfYYrkTDYJ4YC7H2w96GZm/QGVf0qvNvq4AXzq/El86vxJvLZ5NzWpNB9f8CJjhgwMZmGdUcmo4wbkO0yRSGUzZPRPwA7gfwGfAz4DrHH3L+c+vCNpyEjyqbnF+d2bW6lZlubXr2zh7BNGUZ2IcclpYxjQvyTf4Ym0K8pjCP0Irq08j2C20yeBH3se5rxQQpBCsefAIR5ftZGaVJpXN+3m6nAW1lmFWsIqfVrUVUZlwGkEQ0WvuvvBY7wlJ5QQpBCltwclrDWpegCSs2LckIgRHzEoz5GJBKLcQ7gauAt4k2APYRLwZ+7+eBSBdoYSghQyd2f5OzuoTdXz6MoGTh03hGQizlXTxnOcSlglj6JMCK8A8939jfD5icAv8zG/kRKCFIsDh5p5+pUt1KTqeeGtd7nstDEkE3HOO2m0Slilx0WZEF5y9zMznhvBRXLO7OBtOaGEIMXo3T0HeGRFA7V19WzetT8sYY1zikpYpYd0OyGYWTJ8OBeYCPyC4BjCB4G33f0zEcWaNSUEKXavb95NTaqeB+vqGT2kjOSsONfOrGS0Slglh6JICPd09EZ3/2QXY+syJQTpLZpbnN+/+S41qTS/WruZOVUjqZ4d59LTxjCwVCWsEq1Iq4wKhRKC9EZ7DhziiZc3UZtKs2bjLq6aNp7qRJzEBJWwSjSiPIYwieCEtCoyzmyOavrrzlBCkN6ufse+sIQ1TUuLk0zEuWFWjONHqoRVui7KhLACWACsAlpa26OY/rqzlBCkr3B3VqR3UptK88iKBk4ZO4TqRJwrp41jyMDSfIcnRSbKhPCiu58VWWTdoIQgfVFQwtpIbSrN7998l0tOG0MyEeOCkytUwipZiTIhfBQ4GVgEHGhtd/dUd4PsLCUE6eu27T0YlLCm0mzcGZSwJhMxThs3NN+hSQGLMiF8A/gEwZnKrUNG7u6XdjvKTlJCEDnsjS27qU3Vs7CunpGDy0gm4lw7o5KKISphlSNFmRDeAKbka/6iTEoIIu/X3OK88FZQwvrUmqCENZmIc9lklbBKIMoL5LwMDAe2dDsqEYlcST/jvJNGc95Jo9l74BCPv7yJ+5Zs4MsPruLKqeOpTsSYPXGESljlmLJJCMOBV8zsJY48htDjZaci0rHBA/pz4+w4N86O07BjHwvr6rm1ZiWHWpzkrDjJhEpYpX3ZDBld1Fa7yk5FioO7szK9k5pUmkdXbuSkiuOonh3jymnjGaoS1j5BZyqLyPscPNTCM69uoSaV5ndvvMvFrSWsJ42mf0m/fIcnORLlQeXdHL6GchlQCux19x6vc1NCEInO9r0HeXRlAw+k6mnYsY/rZlRSPTvO5PEqYe1tIjuo7O5/mKM3nPr6OuDs7oX3h/VdAXwbKCG4LOedUaxXRI5txOAyPnFOFZ84p4o3G/ewMFXPn967lKHlpVQnYlw3M6YS1j6mU/uIHngQ+EB3OzazEuB7wJXAFOAjZjalu+sV6ZbGRnjppeC+Pb/9Ldx+e3Df1fVk08/atXDvvcF9V2XTT2MjJ65fy18nRrL4by/h7+dPZu3G3Vz6rWf45D1LeGRFA/ubmrvfVxSfJ5t+svnMUfTTG7l7hzcgmXG7EbgT+P2x3pfFes8Bnsx4fhtwW0fvmT17tovkzH33uZeXuw8bFtzfd9/7l5k71x0O3+bN6/x6sunnlluO7OeWW3LzeTpYZu+BJq9NveMf//ELPuOrT/oXa1b4knXvektLS+f7iuLzZNNPNp85in6KDLDUs/lePuYCcE/G7UfAl4Ex2az8GOu9kWCYqPX5J4DvdvQeJQTJmS1bgn/8zC+t8vKgvdXzzx/5euvt+eezX082/axZ03Y/a9ZE+3myWSa0ccc+/8Ezb/jl33rGL/jmb/xfF73q67fuyW49UXyeqLZtVNuuyGSbEI45ZOTun8y4fdrd73D3HjtJzcxuNrOlZra0sS/tuknPWr8eysqObCstDdpbLVrU9nsz24+1nmz6WbKk7X7aa29LNv1ks0xo3LCB/PlFJ7LoCxfyvY8m2LmvieT3f8cH7/od9z+9lp3HDW9/PVF8nmzi7cTn6VY/vVi7B5XN7P928D539691s+964PiM5/Gw7eiO7gbuhqDKqJt9irStqgoOHjU7S1NT0N5q3jz4h394/3vnzct+Pdn0M2dO2zG2196WbPrJZpmjmBnT4sOYFh/Gl66azLOvNVL7+zf5x4/9Oxe9tYzql3/DBetS9M9cTxSfJ5t4u/B5utRPL9bRHsLeNm4ANwG3RtD3S8DJZjbJzMqADwMPR7Bekc6rqIAFC6C8HIYODe4XLAjaW5133pFf/hA8P++87NeTTT+TJ8MttxzZzy23BO1Rfp5slulAWf9+zJ0ylh/cdC6LE82cvfEVvnPBxzjns/fy9a/9N2uaBkT3ebKJt5ufJ+t+erGsTkwzsyHA5wmSwS+Ab0UxbGRmVwH/TlB2+hN3v6Oj5XUeguRcY2MwNFBV1f4XwG9/GwwTHZ0MOrOebPpZuzYYVpkzp/Nfnp3pJ5tlOtHXW0PHUrs+mDZjyMD+VCfiXDezkjH167r/ebKJN+LP0+31FIBITkwzs5HAXwEfA+4Fvu3u2yOLspOUEESKR0uL8+K6bdSk0ixavYnExBEkE3HmTRmrWVh7WLcTgpn9M0Gp6d3A99x9T7Qhdp4Sgkhx2newmUVrNvHAsjQr0zu54vRxVM+Oc2aVZmHtCVEkhBaC2U0PcXjqCgAjOKisqStEpNM27dzPQ8vrqUml2dfUzA2z4lQnYkwcNTjfofVamtxORAqau7O6YRc1qTSPrGigatRgkok4V08fz7ByzcIaJSUEESkaTc0tPPtqI7V1aRa/tpULT6kgmYhx4SkVlGoW1m5TQhCRorTzvSYeWdlAbSrN29ve49oZMZKJGKdXDtXxhi5SQhCRordu615qU2lqU0EJazIR4/qZMcYMHZjv0IqKEoKI9BotLc5L64MS1ide3sTMCSOoTsSYN2Uc5WUqYT0WJQQR6ZVaS1hrUvUsf3s7V0wdRzIRZ07VSPr105BSWyK7QI6ISCEpLyvhupnBBXw27wpKWG9/aDV7DhwimYiRTMSZNFolrF2hPQQRKXqtJawL6+p5aHkDE0aWk0zEuWZ6JcMGqYRVQ0Yi0ic1Nbfw3GuN1NbV89xrjVxw8miSs+JcdGrfLWFVQhCRPm/ne008uqqBhal61m3dyzUzKqlOxJka61slrEoIIiIZ1reWsNbVM6ishGQizvUzY4wb1vtLWJUQRETa0NLiLN2wnZplaZ5YvYnp8WFUJ+LMO30sg8p6Z52NEoKIyDHsb2rmydWbWFhXT2rDduadPo5kIsbZk0b1qhJWlZ2KiBzDwNLDJaxbdu3nweX1/MMja9i9/xA3zIpxQyLGiRXH5TvMHqM9BBGRo6wJZ2F9aHkD8RHlVM+Oc8308QwfVJbv0LpEQ0YiIt10qLmFxa9vpSaV5tnXGjnvxNFUz45zcZGVsCohiIhEaNf+Jh5buZHaVD1vNu7hmhmVJBMxpsWGFXwJqxKCiEiObHh3Lwvr6qlN1VPWvx/JRIwbZsUYP6w836G1SQlBRCTH3IMS1tpUmsdWbWJabBjJRIwrpo4rqBJWJQQRkR60v6mZp9ZspjaVZumG7cybMo7qRIyzT8h/CavKTkVEetDA0hKumVHJNTMq2bJ7Pw8vb+Drv1zLjvcOckM4C2uhl7BqD0FEJIfWbtxFbSrNg8sbiA0vpzoRY/70SkYM7rkS1oIeMjKzDwJfASYDc9w9q295JQQRKVaHmlt4/o2t1KTqeebVLZx74iiSiTiXnDqGsv65LWEt9CGjl4Ek8MM89S8i0qP6l/Tj4lPHcPGpY9i1v4nHV21kweJ13Fa7ivnTx1OdiDM9nt8S1rwkBHdfCxR87a6ISC4MHVjKh86cwIfOnMDb777Hwrp6/uLndZSWBCWs18+MUTm850tYdVBZRCSPJowaxOcvP5m/uOwklm3YTm1dPVd9ZzGnVw6lOhHnA6ePY/CAnvmqztkxBDP7FTCujZe+7O4Phcs8A/x1R8cQzOxm4GaACRMmzN6wYUMOohURKRz7m5r51drNLEzV89L6bVw+ZSxfu25qlxND3o8huPvlEa3nbuBuCA4qR7FOEZFCNrC0hPnTK5k/vZLG3Qf41drNDCoryXm/GjIc6JeQAAAGwUlEQVQSESlgFUMG8JE5E3qkr7xM12dmN5hZGjgH+KWZPZmPOERE5LB8VRktBBbmo28REWlb8UzoLSIiOaWEICIigBKCiIiElBBERARQQhARkZASgoiIAEV2PQQzawS6OnfFaGBrhOHkWjHFW0yxQnHFW0yxQnHFW0yxQvfinejuFcdaqKgSQneY2dJs5vIoFMUUbzHFCsUVbzHFCsUVbzHFCj0Tr4aMREQEUEIQEZFQX0oId+c7gE4qpniLKVYorniLKVYorniLKVbogXj7zDEEERHpWF/aQxARkQ702oRgZuvNbJWZLTezpWHbSDN7ysxeD+9H5DtOaDfWr5hZfdi23MyuynecrcxsuJk9YGavmNlaMzungLdtW7EW5LY1s1MzYlpuZrvM7C8Lcdt2EGuhbtsvmNlqM3vZzO43s4FmNsnMXjSzN8zsf8ysLN9xtmon3p+a2bqMbTsz8n5765CRma0HznD3rRlt/wRsc/c7zeyLwAh3vzVfMWbEtZ73x/oVYI+7/0u+4mqPmd0LLHb3H4f/RIOAL1GY27atWP+SAt22rcysBKgHzgI+SwFu21ZHxfpJCmzbmlkMeB6Y4u77zOwXwGPAVUCtu//czO4CVrj7D/IZK3QY78XAo+7+QK767rV7CO24Drg3fHwvcH0eYylKZjYMuBBYAODuB919BwW4bTuItRhcBrzp7hsowG17lMxYC1V/oNzM+hP8KNgIXAq0frkW2nY9Ot6Gnui0NycEBxaZ2TIzuzlsG+vuG8PHm4Cx+QntfdqKFeAWM1tpZj8phGGC0CSgEbjHzOrM7MdmNpjC3LbtxQqFuW0zfRi4P3xciNs2U2asUGDb1t3rgX8B3iZIBDuBZcAOdz8ULpYGYvmJ8Ehtxevui8KX7wi37b+Z2YCo++7NCeF8d08AVwKfNbMLM1/0YKysUMbL2or1B8CJwEyCP4pv5TG+TP2BBPADd58F7AW+mLlAAW3b9mIt1G0LQDi0dS3w/45+rYC2LdBmrAW3bcOkdB3BD4RKYDBwRV6D6kBb8ZrZx4HbgNOAM4GRQOTDhr02IYRZFnffQnC5zjnAZjMbDxDeb8lfhIe1Fau7b3b3ZndvAX5EEH8hSANpd38xfP4AwZduIW7bNmMt4G3b6kog5e6bw+eFuG1bHRFrgW7by4F17t7o7k1ALXAeMDwckgGIExwHKQRtxXuuu2/0wAHgHnKwbXtlQjCzwWY2pPUxMA94GXgY+ONwsT8GHspPhIe1F2vrF0DoBoL4887dNwHvmNmpYdNlwBoKcNu2F2uhbtsMH+HIIZiC27YZjoi1QLft28DZZjbIzIzDf7NPAzeGyxTSdm0r3rUZPwqM4HhH5Nu2V1YZmdkJBL+0IRg2uM/d7zCzUcAvgAkEs6b+kbtvy1OYQIex/hfBbrcD64E/yxhHzquw3O3HQBnwFkFlST8KbNtCu7F+h8LdtoMJvhBOcPedYVvB/d1Cu7EW5N+tmX0V+BBwCKgD/pTgmMHPCYZf6oCPh7++866deB8HKgADlgN/7u57Iu23NyYEERHpvF45ZCQiIp2nhCAiIoASgoiIhJQQREQEUEIQEZGQEoL0GWbWfNQMnVU57u+nZnbjsZc84j2/C++rzOyjuYlMpG39j72ISK+xz907PWWwmfXPmPMmp9z93PBhFfBR4L6e6FcEtIcgfVz4S3yxmaXC27lh+8Vh+8MEZ7ViZh83syXh3sUPw2mfMbM9ZnaHma0wsxfMLHPyuQvN7Hdm9lbm3oKZ/Y2ZvRROVPbVjPbWE43uBC4I+/pCrreDCCghSN9SnjFc1Hp2+BZgbji54IcIzmJulQA+7+6nmNnk8PXzwr2MZuBj4XKDgRfcfQbwHPDpjHWMB84H5hN8yWNm84CTCeaimQnMPnryRYJJ+Ba7+0x3/7coPrzIsWjISPqStoaMSoHvhlNcNAOnZLy2xN3XhY8vA2YDLwVTyVDO4UnmDgKPho+XAXMz1vFgONHbmow9h3nhrS58fhxBgniuG59NpNuUEKSv+wKwGZhBsMe8P+O1vRmPDbjX3W9rYx1NfngOmGaO/L/KnBvHMu6/4e4/7E7gIlHTkJH0dcOAjeGv+E8AJe0s92vgRjMbA3+4PvfELvb5JPApMzsuXFesdb0ZdgNDurh+kS5RQpC+7vvAH5vZCoKLj+xtayF3XwP8HcGV7VYCTxEcH+i08OpX9wG/N7NVBNdpOPrLfyXQHB6o1kFl6RGa7VRERADtIYiISEgJQUREACUEEREJKSGIiAighCAiIiElBBERAZQQREQkpIQgIiIA/H9fvnyNR3zSMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data (Use Scatter) (does a linear model work here?)\n",
    "df.plot(kind='scatter', x= 'Farenheit', y = 'Number of O-rings damaged', color='r')\n",
    "\n",
    "#Define variables\n",
    "X = df['Farenheit'].values.reshape(-1, 1)\n",
    "y = df['Number of O-rings damaged']\n",
    "\n",
    "# # Fit a linear regression model\n",
    "XX = np.ones((X.shape[0], 2))\n",
    "XX[:,1] = X.ravel()\n",
    "bHat = linalg.inv(XX.T.dot(XX)).dot(XX.T).dot(y)\n",
    "b0, b1 = bHat\n",
    "\n",
    "X_test = np.linspace(50, 85, 100)\n",
    "\n",
    "plt.plot(X_test, b0 + b1 * X_test, linewidth=1);\n",
    "# Conclusion: Nope, does not work!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='multinomial',\n",
       "          n_jobs=None, penalty='l2', random_state=None, solver='lbfgs',\n",
       "          tol=0.0001, verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train the model\n",
    "clf = LogisticRegression(solver='lbfgs', multi_class = 'multinomial') #used for multi-class\n",
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.66928509])"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Prediction: the probability and of P(y = 0 | X = 31)\n",
    "pred = clf.predict(np.array([31]).reshape(-1,1))\n",
    "prob = Sigmoid(pred).ravel()\n",
    "true_prob = (1 - prob) * 100\n",
    "true_prob\n",
    "\n",
    "#The probability of no O-rings fails at 31 degree Farenheit is less than 1%. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://www.davidplopez.com/assets/ChallengerORingData.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
